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ABSTRACT 



Advanced surveillance and weapon guidance systems using 
new mosaic sensor arrays and large scale integration (LSI) 
electronic data processors on the same focal plane require 
detection of weak targets deeply buried in background clutter 
noise by many tens of db ' s . This investigation reports on 
the focal plane processing techniques to accomplish these 
heretofore unachievable goals. 

Five focal plane processing algorithms are developed con- 
sisting of nonrecursive statistical spatial- temporal filters 
for clutter suppression followed by thresholding for initial- 
ization of target detection. These filters are based on either 
the minimization of mean square error criterion (for MMSE 
filters) or the maximization of signal to noise ratio cri- 
terion (for matched filters) . Two are nonadaptive spatial 
filters and two are nonadaptive spatial- temporal filters. 

The fifth type is an adaptive spatial filter based on the 
minimization of mean square error criterion. 

These filters have been investigated analytically and by 
computer simulation using computer generated images contain- 
ing correlated clutter noises modeled by Markov processes and 
also real world infrared images. Using an infrared image in 
the red spike spectral band, a single frame statistical 
spatial filter can suppress the clutter by 27 db. A five 
frame sequential statistical spatial- temporal filter was 
found to have a clutter suppression of 87 db . 
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I . INTRODUCTION 



A. MOSAIC FOCAL PLANE AND SMART SENSOR DEVELOPMENTS 

The enormous impact of integrated circuit technologies on 
baseband electronics is well known. Revolutionary advances 
of LSI (large scale integration) and VLSI (very large scale 
integration) devices are making high performance computer and 
signal processing capabilities available at lower cost with 
lighter weight, less power consumption and higher reliability. 
Their impact on military systems will be tremendous and far- 
reaching . 

In recent years, IC technologies have also been making 
important changes in high frequency electronics throughout 
the microwave and optical bands. In optical detection tech- 
nologies, thousands and hundreds of thousands of optical and 
infrared detectors can be made in the not too distant future 
on one chip by IC batch processing in the form of what is 
generally known as "mosaic detector array" or "mosaic focal 
plane array" [1]. The availability of such a large number 
of detectors on the same focal plane has inspired new modes 
of optical detection, such as the time-delayed- integrat ion 
scanning [2], staring [5], and step-staring [3] modes which 
offer better detection performance. Combined with LSI com- 
puters or signal processors on the same focal plane, complex 
tasks can now be carried out at or near the focal plane by 
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what are called "smart sensors."'*' This equipment is intelligent, 
small, light, consumes little power and is not expensive. It 
is certain that these smart sensors are going to make signif- 
icant improvements to many surveillance, reconnaissance and 
guidance systems. 

This thesis deals with an image processing, or more 
specifically, a "focal plane processing" problem which is 
created by these new "mosaic focal plane" and "smart sensor" 
developments. It will be explained in the next section. 

B. REQUIREMENTS OF NEW FOCAL PLANE PROCESSING 

With the new-found computing and signal processing capa- 
bilities which will be mounted next to the detector array in 
the near future, very ambitious focal plane processing objec- 
tives have been formulated. In front of the mosaic detector 
array, computer controlled adaptive optics and programmable 
optical filters are being developed. On the mosaic detector 
focal plane, sophisticated "on chip" signal processing oper- 
ations such as time-delayed integration [2], a.c. coupling [4], 
differentiation [5], multiplexing [6], convolutional read-out 
and others are being investigated. After, or during, the 
read-out, image processing operations of sampled analog signals 
can be carried out. If digital focal plane processing is used, 
signal formating, sample and hold and analog to digital 
conversion must he performed. 

■*■ "Smart sensors" session in SPIE 1979 East Symposium to 
be held in Washington, D.C. April 17-18, 1979. 
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For military applications of surveillance, reconnaissance 
and guidance, focal plane processing objectives include the 
basic signal manipulations described above and other predetec- 
tion steps, such as background clutter suppression, target 
enhancement and thresholding for initialization of target 
detection. Following this, there are many post - thresho lding 
focal plane processing steps such as target linking, track- 
ing, track assembly, track association, trajectory estimation, 
target discrimination, target recognition, and many others. 

This thesis is concerned with the predetection steps of 
background clutter suppression and thresholding for initial- 
ization of target detection. 

In visible imaging systems, the contrast between targets 
and background is usually high. The target signal to noise 
ratio is often quite adequate. Therefore, focal plane process- 
ing operations of visible images can directly deal with target 
signals and background noise in most cases without any pre- 
processing to improve the signal to noise ratio. However, 
in infrared imaging systems, the contrasts between targets 
and backgrounds are much smaller. Target signals are often 
a very small fraction of the incoming signal which consists 
of a large amount of background clutter noise and other white 
noise. All infrared imaging systems must deal with the 
clutter problem. The need of suppressing background clutter 
and detecting targets buried in clutter certainly is not unique 
in these new mosaic focal plane arrays. However, the 



11 



requirements of clutter suppression and target detection are 
now much more demanding because of the improved performance 
of mosaic focal plane sensors processors and much more 
ambitious system goals. In new surveillance and guidance 
systems, it is sometimes required to detect weak targets 
deeply buried in clutter by many ten's of db ' s . Old pro- 
cessing techniques used to separate targets from clutter are 
simply not adequate to meet these requirements. New process- 
ing techniques using spatial, temporal and spectral character- 
istics of both the targets and clutter must be developed. 

These techniques are different from most of the well known 
image processing techniques in two major aspects: 

° They deal with image signals of very low signal to 
noise ratio. Furthermore, the dominant noise com- 
ponent is due to clutter and is correlated rather 
than white. 

0 The signals of interest are usually point targets, 
line targets and targets of different shapes instead 
of pictorial details of the image. 

Therefore, they are called "focal plane processing" in this 
thesis instead of the widely used name of "image processing" 
which usually concerns pictures, photographs, TV pictures, 
etc . 

C. FOCAL PLANE PROCESSING TECHNIQUES FOR BACKGROUND 

CLUTTER SUPPRESSION AND TARGET DETECTION 

Clutter suppression and target detection are not unique 
problems for mosaic focal plane processing. Similar types 
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of problems exist in radar and sonar processing and have been 
extensively investigated [7,8]. Effective processing algor- 
ithms for clutter suppression and thresholding for target 
detection have been well developed. However, there are sig- 
nificant differences between mosaic focal plane systems and 
radar and sonar systems so that most of the developed process- 
ing algorithms either are not applicable or must be modified 
if they are to be used for imaging systems. The major dif- 
ferences are the following: 

1. Imaging systems are mostly passive. Therefore, 
processing algorithms based on Doppler shift or the signal 
waveform cannot be applied to passive optical systems. 

2. Image signals are three dimensional, which includes 
two dimensional spatial variations on a single image frame 
and one dimensional temporal variations on several frames of 
images. Radar and sonar signals are 
with temporal variations. Therefore, 
cessing algorithms must be modified 1 
to mosaic focal plane problems. 

In spite of these differences, th 
processing approaches can he applied 
the following, several processing approaches for target 
enhancement and detection are listed. 

When signal to noise ratio is above unity, thresholding 
can be applied to separate true targets from background noise. 
The threshold level can be determined from several different 
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considerations; for example, based on a constant false alarm 
rate (CFAR) , or on a given probability of detection, or on a 
given probability of tracking, etc. 

If the signal to noise ratio is low, as in many practical 
situations, prethresholding signal processing steps must be 
employed to suppress the noise and enhance the target, if 
possible, so that the signal to noise ratio will be adequate 
for thresholding. Some of these processing approaches are 
based on the following concepts: 

Summing - Time integration 
Differencing - pulse canceller 
Filtering : 

Based on time domain knowledge - correlator 

Based on frequency domain knowledge - frequency 
domain filter 

Other approaches are statistical, based on the following 
procedures : 

Statistical processing: 

Maximization of signal to noise ratio - statistical 
matched filter 

Minimization of mean square error 
Kalman filter 

Wiener filter, MMSE filter 
Maximization of likelihood filter 

The principles of some of these signal processing approaches 
have been extended and applied to the focal plane processing 
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background clutter sup- 
rocessing, based on the 

Focal Plane Processing: 

Spatial filters : 

Second order (Laplacian) 
[9] 

integration [2] 

- Temporal filters: 



9 



Summing - Multiple frame integration 

Differencing - Frame to frame differencing. Non- 
recursive Second differencing, Third 
differencing, etc. [10] 

Filtering - Recursive digital high pass filters 
based on the knowledge of frequency 
characteristics of the clutter [11] 
Active band pass filters. 



Several statistically designed processing approaches have 
also been investigated. 

Statistically Designed Focal Plane Processing: 

Single frame processing: 

Spatial Kalman filter [12,15] 

Spatial MMSE filter [14] 

Multiple frame processing - Temporal filter 
Maximum likelihood filter [15] 

Although quite a few processing algorithms are listed above, 
only the deterministic single frame spatial filters and the 
multiple frame differencing filters have been substantially 
developed. All other processing approaches have only received 
preliminary investigation. None of these filters have shown 
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enough clutter suppression capabilities to enable detection 
of weak targets more than 70 db below clutter. 

D. FOCAL PLANE PROCESSING TECHNIQUES 

DEVELOPED IN THIS THESIS 

The objective of this thesis investigation was to develop 
statistical processing algorithms for background clutter 
suppression. This is accomplished in several phases. First, 
design procedures are developed for two types of nonrecursive 
spatial filters. These concepts are then extended to three 
dimensions by adding a temporal filter. Several aspects of 
target detection are also addressed, with special considera- 
tion given to cases involving non-Gaussian probability 
densities. Adaptive filter algorithms are also developed 
for situations involving unknown and non- s tationary statistics 
of background clutter noise. Analytical results and computer 
simulation studies demonstrate the feasibility of the filter 
algorithms proposed. The following paragraphs describe in 
more detail the development of the processing techniques. 

In Chapter II the basic structure of nonrecursive spatial 
and/or temporal filters is introduced. The spatial filters, 
which may be adaptive, utilize image signals in the neighbor- 
hood of the "estimation pixel" to suppress clutter and enhance 
target intensity at the estimation pixel. Further improvement 
is obtained by also processing spatially filtered data from 
adjacent image frames. 
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Expressions for the filter coefficients are derived in 
Chapter III for two criteria which have been previously used 
in one-dimensional statistical processing filters [16]: min- 

imization of mean square error (MMSE) and maximization of 
output signal to noise ratio (matched filter). Statistical 
nonrecursive spatial filters based on the MMSE criterion have 
been studied before for image restoration and for improvement 
of signal to noise ratio against white noise [23-25]. The 
design of these filters, however, is carried out using spatial 
frequency domain techniques, or equivalent. Such filters have 
not been very useful in practical applications because of the 
large size, e.g., larger than 10x10, required to achieve 
meaningful improvement [25]. Furthermore, a nonrecursive 
spatial filter tends to blur the edges of an object [17]. The 
concept of using nonrecursive statistical spatial filters for 
suppression of correlated noise was originated by Bar Yehoshua 
[14] who developed a design procedure in the spatial domain 
using the MMSE criterion. Bar Yehoshua developed the special 
application of statistical nonrecursive spatial filters for 
enhancement of point/line target signal to noise ratio in the 
presence of spatially correlated noise. He showed that spa- 
tial filters as small as 3 x 5 pixels are effective in the 
suppression of spatially correlated noise. Furthermore, his 
design procedure in the spatial domain is well suited for hard- 
ware implementation because the rather complex two dimensional 
transformation operations required in spatial frequency 
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domain design procedures are unnecessary. Chapter III reviews 
Bar Yehoshua's work and extends the concept to nonrecursive 
spatial matched filters. Other related work is contained in 
[17] . The matched filter for the enhancement of point and 
line targets described here is applicable to correlated, as 
well as white, noise. In addition, a new explicit design 
procedure is presented. This thesis extends Bar Yehoshua's 
concept in three directions. In this chapter, a different 
design criterion based on the maximization of signal to noise 
ratio is developed for the second type of nonrecursive spatial 
filter- -matched filter. Again, spatial matched filters are 
not new [17,13]. The contribution of this study is in the 
development of its new application to the enhancement of point/ 
line target to noise ratio in two dimensional and three dimen- 
sional images of very low signal to noise ratio, and in the 
design procedure in spatial and temporal domains instead of 
the spatial frequency and temporal frequency used by most 
other workers. Extensions in two other directions will be 
described in Chapter IV and Chapter VII. 

In Chapter IV, two nonadaptive multiple frame processing 
approaches for clutter suppression are developed. They are 
the nonrecursive sequential spatial - temporal filters designed 
again by two statistical criteria similar to those in Chapter 
III which lead to the spatial -temporal MMSE filters and the 
spatial- temporal matched filters. Temporal filters which 
use several frames of image designed by either deterministic 
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or statistical information are not new [17]. However, the 
combinations of spatial and temporal filters using statistical 
information of several frames of image for achieving sub- 
stantial enhancement of signal to noise ratio has not been 
reported before. This thesis investigates a special type of 
multiple frame processing consisting of sequential applica- 
tion of a spatial filter followed by a temporal filter using 
several frames. It is shown analytically that the clutter 
suppression of the spatial filter and of the temporal filter 
are multiplicative (or additive, if they are expressed in 
terms of db) . 

In Chapter V, the thresholding problem for initialisation 
of target detection is investigated. The concept of thresh- 
olding for detection has been well developed for radar 
systems [8]. The extension of the thresholding operation to 
two dimensional problems has not been extensively investigated for 
cases where the probability density functions are not Gaussian. 
In this chapter, a procedure for setting the threshold based 
on the assumption of a Gaussian distribution is applied to 
computer generated images whose probability density functions 
are Gaussian. In addition, a new procedure for setting 
threshold level based on a constant false alarm rate criterion 
without the assumption of a Gaussian probability density 
function is also developed. 
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In Chapter VI, the four nonadaptive statistical filters 
developed for clutter suppression are combined with thresh- 
olding as the nonadaptive focal plane processing algorithm 
studied in this thesis. Their properties and performance 
are quantitatively studied by computer simulation using two 
types of images. The first type consists of computer gener- 
ated images which contain targets and two random noise com- 
ponents. One noise component is uncorrelated, or "white" 
Gaussian noise. The other noise component is correlated, or 
"colored" noise modeled by two Markov processes- - first order 
and second order [14,19]. These correlated noises are used 
to simulate background clutter. In these computer generated 
images, the statistical properties of the correlated noise 
are known "a priori." The second type of image is a set of 
twenty frames of real world infrared images. Here, correlated 
noise in both spatial and temporal domains is not known "a 
priori." Practical filter design procedures are developed 
which first "learn" the statistics from the real world image, 
then design a filter using the techniques developed in Chapters 
III and IV. Using single frames of these real world images, 
clutter suppression of up to 27 db has been achieved by spatial 
filtering alone. Using five frames, clutter suppression of 
up to 87 db has been achieved by sequential spatial- temporal 
filters. These large clutter suppression values based on 
computer simulated tests on real world infrared images are 
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the most important testimony of the practical usefulness of 
the statistical focal plane processing algorithms developed 
in this thesis. 

In Chapters VII and VIII, adaptive design procedures are 
developed to further extend the concept of clutter suppression 
using the nonrecursive spatial filter. No "a priori" knowledge 
of the statistical properties of the clutter is needed to 
design the filter. The adaptive procedure automatically 
corrects the filter coefficients iteratively toward the 
optimal values to yield the minimum error between the target 
estimate and a selected desired response. The concept of one 
dimensional adaptive filters based on the least mean square 
error criterion has been very well developed by Widrow, McCool, 
et al. [20-22]. This thesis extends the basic concept to two 
dimensions and also develops the following new techniques to 
improve the performance of the adaptive spatial filter: 

1. Derivation of desired responses corresponding to 
various types of targets of interest. 

2. An "automatic gain control" procedure which yields 
constant misadj us tment for varying statistical properties of 
the clutter, prevents overflow and assures adaptive conver- 
gence . 

3. Dividing an image into small subimage tones to carry 
out filter adaptations if the statistical properties of the 
image are both "nonstationary" and unknown a priori. 
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Simulation results are given in both Chapters VII and 
VIII which demonstrate the efficacy of the proposed design 
procedures . 

Chapter IX contains a summary of the results developed 
in this thesis and conclusions. 



II. BASIC CONCEPT AND PROCEDURES OF NONRECURSIVE 
STATISTICAL FOCAL PLANE PROCESSING 



A. INTRODUCTION 

Focal plane processing has a wide range o£ objectives 
such as optical filtering, sensing, multiplexing, signal 
formating, background clutter suppression, thresholding, 
target tracking, target discrimination, target recognition, 
and many others. This thesis is only concerned with back- 
ground clutter suppression and thresholding. They are two 
old problems taking on new urgency and demanding requirements 
due to the recent developments of "mosaic focal plane arrays" 
and "smart sensors." 

The thresholding process is not significantly different 
from that developed for radar. However, clutter suppression 
algorithms will be drastically different from those developed 
for radar and for old infrared systems using only a small 
number of detectors. New focal plane processing techniques 
capable of suppressing clutter by several ten's of db ' s are 
urgently needed to meet demanding requirements of new sur- 
veillance and guidance systems. 

Five nonrecursive statistical focal plane processing 
algorithms are developed in this thesis for clutter suppression. 
These are followed by a thresholding step to initiate the 
target detection process. Although differing in detailed 
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implementation, they share the same basic concept and filter 
configurations. The main difference lies in the design of 
their filter coefficients. They will be presented in this 
chapter. An appreciation of the common underlying concept 
will help to place in pr oper • per spective the detailed devel- 
opments of all five statistical filters in the following 
chapters . 

B. BASIC CONCEPT OF NONRECURSIVE STATISTICAL 

FILTERS FOR CLUTTER SUPPRESSION 

One of the major objectives of focal plane processing 
developed in this thesis is background clutter suppression. 

In infrared mosaic focal plane array systems, clutter origin- 
ates from thermal radiations in the background scene in the 
field of view of the sensor system. Since the background 
scene has well defined spatial shapes and temperature distri- 
butions, the background clutter noise in the image must have 
well defined spatial shapes or structures also which may be 
contaminated by scattering through the atmosphere, distortion 
through the optical and electronic components, and different 
noise sources. In other words, the clutter noise is a random 
signal with well defined spatial structure for a given scene. 
This information on the clutter noise is useful in helping to 
achieve its suppression. However, the important concept of 
the statistical filter is the assertion that although the 
spatial structure of the background clutter varies with the 
scene, the statistical characteristics of the spatial structure 
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over a large area may be described by well defined statis- 
tical characterizations such as the correlation matrix, 
autocorrelation function or power spectral density. If this 
is true, a spatial filter can be designed to use the signals 
around a selected pixel or over its neighborhood to process 
the signal at that pixel. Since background clutter has well 
defined statistical spatial behavior while the target of 
interest only dwells in one pixel, or several pixels of well 
defined location, signals in the neighborhood can be properly 
used to suppress the spatially related background clutter and 
thus raise the target signal out of the clutter. This is the 
basic concept of the statistical spatial filter for clutter 
suppression. The filter can be designed by performing a 
statistical estimation of the target signal, using one of 
several criteria. They will be explained in a later section. 

The statistical correlation concept in the spatial domain 
can be extended to the temporal domain also. If successive 
frames of images are fairly well registered, with small amounts 
of drift and jitter in between frames, it is reasonable to 
expect that the signals are correlated in temporal domain. 

If this is true, a temporal filter can be designed that will 
use the signals at the same pixel location of successive frames 
of images to process the signal at one of the frames for sup- 
pression of the temporally correlated clutter. 

In summary, if background clutter sources are related 
(in the statistical sense) either spatially in a single frame 
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or temporally over several frames, signals in both the spatial 
and temporal neighborhood can be used together to process the 
signal at the selected pixel for suppressing the clutter. 

These unprocessed signals are multiplied by the filter coef- 
ficients and summed together to provide the processed signal 
at the selected pixel. Therefore, it is a "nonrecursive" 
filter operation. The filter configuration will be described 
in the following section. 

C. CONFIGURATION OF NONRECURSIVE FILTERS 
1. Nonadaptive Nonrecursive Filters 

In the previous section, the concept of nonrecursive 
filtering using unprocessed signals in the neighborhood around 
a selected pixel, to be called the "estimation pixel," was 
introduced. The purpose of such a filter is to suppress the 
background clutter by performing a statistical estimation of 
the target signal at that pixel. The neighborhood area where 
the filtering is carried out is called the "search box." The 
location of the "estimation pixel" is usually inside the 
"search box" but can also be outside. For symmetrical clutter 
characteristics, it is best to have the "estimation pixel" in 
the middle of the search box. However, it must be moved 
toward the edges of the search box when the estimation pixel 
is located at the edge of an image because there is no signal 
outside the image. 
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In summary, unprocessed signals in a search box are 
used to process the signal at only one "estimation pixel" 
inside the search box. The estimation pixel is then scanned 
throughout the image. This filter configuration is further 
described in two cases shown in Fig. 2-1 and 2-2 for a spatial 
filter and a spatial - temporal filter respectively. 

T D T T 

Let X = unprocessed signal = S + N (2-1) 

T D 

S = targets of interest = [s-, s 7 s. ...s,](2-2) 

a , - , l , l 

t D 

N = sum of all noise components including (2-5) 
uncorrelates "white" noise and 
correlated "colored" noise = 

[ n l , n 2 , n i , * * - n i^ 

where i = 1,2,... I and I is the total number 

of pixels in the search box. 

CASE I - Spatial Filter (Fig. 2-1) 

A nonrecursive spatial filter is described by a set of 
filter coefficients, or weighting coefficients, over the area 
of a "search box." The set of filter coefficients is repre- 
sented by a vector W. The intensity of the unprocessed 
signals in the search box, represented by the vector X, are 
multiplied by the weighting coefficients, (the vector 'Tv 1 ') corre- 
sponding to the same pixel location in the same search box and 
summed together to give the output intensity at the estimation 
pixel, y, that is, 

y = W T -X 

The filter coefficients are designed by different criteria to 
be described in the following section. Spatial-filtering using 
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SPATIAL FILTER 




Figure 2-1 



Configuration of Nonadaptive Two Dimensional 
Tonrecur s ive Spatial Filter 



one search, box will yield a processed signal, called the 
"estimation signal,." of only one pixel, the "estimation pixel." 

The search box is then moved, following the scanning of the 
"estimation pixel" throughout the image. 

This thesis concentrates on a search box size of 3 by 5 
pixels. Therefore, nine filter coefficients are required. 

The estimation pixel is usually the center pixel. It is moved 
to the middle of one edge or one corner sometimes, as explained 
above . 

CASE II - Spatial -Temporal Filter (Fig. 2-2) 

In general, a three dimensional nonrecursive filter is 
described by a set of filter coefficients, W, over the volume 
of a "search box" which consists of several frames of a two 
dimensional box. For example, in a 3x3x5 search box, a 3x3 
search box is used in five successive frames of the image. 
Altogether, W consists of 45 coefficients. However, the 
spatial- temporal filters developed in this thesis are the 
special case which consists of sequentially applying a spatial 
filter, first on one frame of image, then followed by a tempo- 
ral filter, W^, using several successive frames of images. 

For example, consider a sequential 5x5 spatial filter followed 
by a five frame temporal filter. The spatial filter is imple- 
mented by the same procedure described in Case I. For the 
five frame temporal filter, signals at the same pixel location 
from five successive frames are multiplied by the temporal 
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Three Frame Sequential Spatial-Temporal Filter 
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Figure 2-2 Configuration of Nonadaptive Three Dimensional 
Nonrecursive Sequential Spatial-Temporal Filter 
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filter coefficients and summed together. Therefore, instead 
of 45 filter coefficients for the general three dimensional 
filters, the sequential spatial- temporal filters have only 
14 filter coefficients. The first nine coefficients are for 
the spatial filter. The last five coefficients are for the 
five frame temporal filter. 

2 . Adaptive Nonrecursive Filters (Fig. 2-5) 

In this thesis, only an adaptive spatial filter is 
developed. Its configuration is an extension of the spatial 
filter configuration shown in Fig. 2-1. The concept of adap- 
tation is based on the one dimensional adaptive filter orig- 
inally developed by B. Widrow [20,21,22]. 

In the adaptive filter, there is an automatic correc- 
tion procedure for adapting the filter coefficients by using 
a feedback loop shown in Fig. 2-3. The search box and oper- 
ation of an adaptive spatial filter are still the same as that 
of the nonadaptive spatial filter. However, the output of 
the spatial filter is now compared with a desired reference 
signal, dj , by a subtraction operation. The error of this 
comparison is multiplied by an adaptive loop gain, y , and the 
unprocessed signal at a pixel to yield a correction term for 
the filter coefficient corresponding to that pixel. This 
adaptive correction is carried out by scanning the search 
box over an "adaptation area." Therefore, in the adaptive 
filter case there are three image areas: 

The given image area. 
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Adaptive Spatial Filter 



Search 

Box 




Unprocessed 

Image 
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Figure 2-5 



The spatial filter area - search box area, 

The adaptation area - which is in general 
larger than the search box area; for example, 
8x8 for a search box of 3x3. 

D. PRINCIPLE OF STATISTICAL FILTER DESIGN 

In section C, configurations of all five statistical 

filters are shown in three different types: Figures 2-1, 

2-2 and 2-3. The common principle behind them is to "sum" 

the "weighted" signals in the neighborhood of the "estimat 

pixel." The neighborhood is defined in a general sense to 

include both spatial neighborhood and temporal neighborhoo 

Once the configuration is chosen, the only parameters rema 

ing to describe the filter are the weighting coefficients. 

They are designed by -several different criteria. 

Two criteria are used in this thesis. 

1 . Criterion I - Minimization of Mean Square Error 

In this criterion, a desired result is selected. 

example, in the nonadaptive spatial filter case, the objec 

of spatial filtering is to estimate the target, s. However 

the output of the filter only provides an estimate of the 

target, y = s. Therefore, there is an error, e = s - s. 

2 

mean square value is E(e“). The filter coefficients, W, a 
designed by minimization of the mean squared error. This 
filter is called "spatial MMSE" filter or "spatial Wiener 
filter" [14,18,23,24,25,26] . 

A similar procedure is developed for the temporal 
MMSE filter and the spatial- temporal MMSE filter. 
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2. Criterion II - Maximization of Signal to Noise Ratio 

[27,28,29] 

In this criterion, the objective is to maximize the 
output signal to noise ratio. The filter designed by this 
criterion is called the "spatial matched filter." 

A similar procedure is developed for the temporal 
matched filter and the spat ial - temporal matched filter. 

E. IMAGES USED IN ANALYZING AND TESTING OF THE 
STATISTICAL FOCAL PLANE PROCESSING 

Two dimensional infrared images containing a background 
scene and different targets have been collected in the past 
using infrared scanners and radiometers. Although they 
provide some information on the background and target, they 
are found inadequate because these instruments typically use o 
a few detectors. They are scanning instruments with relatively 
large instantaneous fields of view and fairly broad spectral 
bands. Many of them are not in spectral bands suitable for 
target detection in the new surveillance systems. Therefore, 
there is an acute lack of standard background scenes needed 
for developing and testing new focal plane processing algor- 
ithms to be used in these new high performance surveillance 
and guidance systems using smart sensors. 

In this thesis, a set of real world infrared images of 

twenty frames with relative drift of 10% pixel each frame, 

? 

selected by a responsible agency, was made available. - 
? 

M. Schlesinger, Aerospace Corporation., private 
communication. 



They are used extensively in testing both the spatial filters 
and spatial- temporal filters. The performance results are 
the most convincing evidence of the usefulness of this research 
because the test is based on real world infrared images. 

However, to carry out detailed research of this statis- 
tical focal plane processing, a wide variety of images with 
varying parameters is needed. Since there is simply no such 
supply of real world images taken by the new mosaic focal 
planes, it was necessary to develop a theoretical model for 
the background clutter. If it can be proven valid, different 
test images can be generated on a computer readily by simply 
varying the appropriate parameters. 

For this purpose, contaminated images have been generated 
by using the following models. A contaminated image consists 
of three components. 

1 . Target 

It is described by its spatial shape, its movement 
between frames and its intensity. The intensity is described 

either deterministically by S, as defined by eq. 2-2, 

? 

or statistically by a mean and variance 

2 . Random Noise 

There are two types. The first type is Gaussian un- 
correlated, or "white" noise described by its mean y^ and 

9 

variance a^. - . The second type is Gaussian correlated, or 
"colored" noise - modeled by two Markov processes [Appendix A]. 
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A first order Markov process described by its auto- 
correlation function. 



_ 2 m 

R " a CN P H 

where = 

°H = 

p V = 

m = 

n = 




variance 

horizontal correlation factor 
vertical correlation factor 
horizontal coordinate of image 
vertical coordinate of image. 



A second order Markov process described by its 
correlation function. 
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horizontal spatial frequency 
vertical spatial frequency. 
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Ill . SINGLE FRAME FOCAL PLANE PROCESSING FOR 

CLUTTER SUPPRESSION - NONRECURSIVE STATISTICAL 
SPATIAL FILTERS 

A. INTRODUCTION 

In this chapter, two statistical spatial filters which 
use only one frame of an image are developed. One is designed 
by minimizing a mean square error criterion. This type of 
filter is called "MMSE" filter, which is more commonly known 
as a Wiener filter. It provides an "optimal" estimate of the 
target signal and will be developed in Section C. The second 
type is designed by maximizing the output signal to noise 
ratio criterion and is commonly known as the matched filter. 

It will be developed in Section D. 

These two statistical design procedures can be equally 
well applied to focal plane processing using several frames 
of images in the form of three dimensional spatial- temporal 
or temporal - spatial filters. They will be developed in the 
next chapter. 

Both the spatial filters developed in this chapter and 
the spatial - temporal filters developed in the next chapter 
will be evaluated by using two types of images. One type is 
computer generated images containing correlated noise modeled 
by a Markov process to simulate background clutter. The other 
type is a real world infrared image. The results from both 
simulated and real world scenes are presented together in 
Chapter VI. 
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B. CONFIGURATION OF TWO DIMENSIONAL 
NONRECURSIVE SPATIAL FILTERS 



Although the design criteria of two statistical spatial 
filters developed in this chapter are different, they share 
the same filter configurations as shown in Fig. 3-1. There- 
fore, the filter configuration and operation procedure will 
be described before their individual design criteria are 
presented in the next two sections. 

In these spatial filters, signals in the neighborhood 
of a selected pixel of interest are used together to help the 
processing of signals at that pixel. The selected neighbor- 
hood is called a "search box." As explained in Chapter II 
and the previous section, the objective is to suppress the 
background clutter by performing an estimation of the target 
signal at that pixel, using its neighborhing signals and the 
knowledge of their statistical spatial correlation properties. 
This is carried out by designing a set of weighting coefficients, 
or filter coefficients, W. For every pixel in the search box, 
a weighting coefficient is assigned. Signals of all pixels, 

X, in the search box are multiplied by their corresponding 
filter coefficients, W, and summed together to give an "optimal" 
estimate of the target intensity at the "estimation pixel," 
which can be anywhere inside the search box. It could also 
be located outside the search box although this is not devel- 
oped in this thesis. In general, the location of the "esti- 
mation pixel" is in, or close to, the middle of the search box 
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SPATIAL FILTER 




Figure 3-1 Configuration of Nonadaptive Two Dimensional 
Nonrecursive Spatial Filter 
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because common sense suggested that spatial behaviors in all 
directions around the estimation pixel should be used to 
reveal the spatial behavior of background clutter and to help 
achieve their suppression. However, there are situations 
where signals in certain parts of the neighborhood do not 
exist; for example, when the estimation pixel is located at 
the edge of an image. In these cases, it is not possible to 
have the estimation pixel in the middle of the filter. 

Instead, it must be off to the side, as will be shown later 
by several examples. 

The size and shape of the search box can also vary. Only 
square shaped search boxes are developed in this thesis be- 
cause the correlation properties in vertical and horizontal 
directions are assumed to be quite close. On the other hand, 
the concept and design procedure developed in this thesis can 
be easily extended to non- square- shaped search boxes if needed; 
for example, if the correlation properties in the vertical 
direction are signif icantly different from those in the 
horizontal direction. 

The size of the search box depends on the correlation 
property of the background clutter and on the size of the 
targets of interest. If the spatial correlation of the clutter 
is strong, there is no need to go beyond a few neighboring 
pixels to determine the correlation properties. Therefore, 
a small search box is sufficient to accomplish an effective 
suppression of background clutter. On the other hand, if 
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the spatial correlation of clutter is weak, one needs further 
away neighboring pixels to obtain the correlation properties 
properly and a larger search box is needed. In this thesis, 
a 3x3 search- box is developed almost exclusively because it 
was found that the correlation properties of several real 
world infrared images are quite strong. In terms of a corre- 
lation factor, on a scale from 0 to 1 corresponding to no 
correlation (white noise) to perfect correlation (determin- 
istic), the real world images examined in this thesis indicated 
a range of values above 0.8. On the other hand, the concept 
and design procedure developed in this thesis can be easily 
extended to large sizes such as 4x4, 5x5, etc. 

In summary, nonrecursive spatial filters developed in 
this chapter consist of a search box for every pixel of in- 
terest in the image, called the "estimation pixel." Because 
secondary signals in the search box are used to estimate the 
target at the estimation pixel, target estimation pixel 
can be moved to any part of the image or scanned to cover any 
selected corner of the image. In this thesis, a raster scan 
of the estimation pixel is used to cover the whole image. 

When the estimation pixel is at the corner or along the edges 
of the image, it cannot be located in the middle of the search 
box. On the other hand, whenever possible it will be located 
in the center. For every search box, signals of all pixels 
are weighted and summed to give an estimate of the target 
intensity at only one pixel inside the box, the estimation pixel. 
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The weighting coefficients are designed using different 
criteria. Two of these design procedures are presented in 
the next two sections. 



C. NONRECURSIVE SPATIAL MMSE FILTER 

1 . Principle of Optimal Estimation 

The objective of the spatial MMSE filter is to give 

a processed signal Y = fy(m o ,n Q ) which will be the "optimal" 

estimate of the target of interest, S = fg(m ,n ) where 

(m ,n ) are the indices of the estimation pixel. The cri- 
v o o J r 

terion for "optimality" is the "minimum mean square error" 
(MMSE), sometimes known as the Wiener filter. The theoretical 



derivation of the MMSE filter is well known and more detailed 
information can be found in the following references: 

- Reference IS states the general MMSE criteria. 

- Reference 23 describes an image restoration technique 
using spatial frequency domain design technique. 

- Reference 24 describes an image enhancement technique 
based on the MMSE algorithm using, again, spatial fre- 
quency domain design method, requiring efficient trans- 
formation techniques. 

- Reference 25 describes in detail the spatial domain 
derivation algorithm for the MMSE filter considering 
some implementation methods in the frequency domain. 

- Reference 26 describes the linear algebraic image restor- 
ation method in general, and the least squares technique 
including the Wiener filter formulation in particular. 

The design equations, both in spatial and spatial fre- 
quency domain, are presented. 
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A practical method of implementation of spatial 
Wiener filter for suppression of correlated noise and enhance- 
ment of point and line target was suggested in Ref. 14. This 
design method will be further extended here and in Chapter IV, 
which presents multiple plane spatial temporal filtering. 

The design procedure is presented in the following paragraphs. 



[Step I] Description of Unprocessed Image: 

The unprocessed signal of the whole image X is given 
by a two dimensional array shown in Fig. 5- 2 (a). 



where 
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intensity of 


unprocessed signal 


(3-3. 1(a) ) 


S = 


target signal 


of interest 




(3-5. 1(b)) 



N = sum of all noise sources (3-3. 1(c)) 

MMAX - number of horizontal pixels in the image 
NMAX - number of vertical pixels in the image. 



The original signals inside the spatial filter search box are 
given by a smaller two dimensional array. 
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where 



m = horizontal pixel number of the upper left 
corner of the search box 

n = vertical pixel number of the upper left 
corner of the search box 
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( 3 - 3 . 3) 



MSB • NSB^ 

This two dimensional array is further simplified by using a 
one dimensional vector representation. 

• • x T ] (3-5.4) 



V T D 
X = [x 



X. 



. X . 



where 



I ’ i * ^ T 

i = MSB • (An- 1) + Am 

A m = Horizontal pixel location in search box with 



n 



respect to top row 

= Vertical pixel location in search box with 
respect to left-most column 



I = MSB MSB - total number of pixels in the 
search box 

The objective of the spatial MMSE filter is to give the 
optimal estimate of the target signal of interest at the esti- 
mation pixel (MES,NES) by summing all weighted signals in the 
search box. (MES,NES) is the location in the search box coordi- 
nate system. In the one dimensional vector representation, the 
weighting coefficients are also described by a vector W as shown in Fig. 5-2 (c). 



4 5 



The Filter Configuration 
Input Samples Vector 




L 



ns u 



C3 -H 



r> 



C/3 U. 



46 



Figure 3- 2(c) Operations of Nonrecursive Spatial Filter. 



W ES, * * * W I^ 



(3-3.5) 




w, w_ ... 

1 , z , 



where w pq = filter coefficient at the estimation pixel 

(mes,nes). 

(Step II] Filtered Image and Mean Square Error: 

The output of the filter is: 



y = X^* W = s (5-3.6) 

/N. 

where s = estimate of the target signal at the estimation 
pixel (MES , NES) . 

Let us define the estimation error as 

e (MES, NES) = s (MES, NES) - s (MES, NES) (3-3.7) 

The mean square error is then the expectation value 

E(e 2 (MES, NES) ] = E{ [s (MES ,NES) - s (MES, NES)]} 

= E { [ s (MES , NES) - X T *W] 2 } (3-5.8) 



[Step III] Minimization of Mean Square Errors: 

To minimize (3-5.8) , it is differentiated with respect 
to W and is then equated to zero. 

JL E ( e 2 (MES , NES) ] = E{ -X- 2 [s (MES ,NES) - X T • W] } 

= E{-X- 2 e } = 0 (5-3 . 9) 

Equation (3-3.9) shows that the estimation error is orthogonal 
to the unprocessed signal X. Further simplification of (3-3.9) 
g ives : 

E (X • s (MES , NES) ] =h[X-X T -W] = E[X«X T ]*W (3-3.10) 

which is the equation for designing the filter coefficient W. 
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(5-3.11) 



Let us define 
D 

P = E[X. s(MES,NES)] 



= cross correlation between the unprocessed 
image signal X and the target signal of 
interes t . 



R xx ^ E[X. X r ] = 



E(x^ x x ) E(x^ x 2 ) . . . E(x-^ Xj) 
E(x2 x-^) E(x2 X 2 ) • • • E( x 2 x I^ 



E (x 1 x ± ) 



E(x i , X P 



(3-3.12) 



= correlation matrix of the unprocessed 
image signal. 



[Step IV] Equation of Filter Design: 

Equation (3-3. 10) is now expressed as: 



R xx -W = P 



or 




(3 - 3 . 13) 



It is a set of I simultaneous linear equations used 
to solve for the weighting coefficients, or filter coeffi- 
cients, of the two dimensional spatial MMSE filter for target 
estimation. 

In many practical applications, we assume that the 
target signal S and the noise N are uncorrelated. If we 
further assume that the mean of the noise is zero, then 
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R xx R ss + r nn 



(5-3.14) 



where 



Rgg = correlation matrix of the target signal 



R 



NN 



= covariance matrix of the noise. 



The filter coefficients (3-3.13) a te given by: 




(3-3.15) 



The subscript WF stands for ’’Wiener Filter." 



2. Design Procedure for Spatial MMSE Filter [Ref. 1] 

To design a spatial MMSE filter, a three step procedure 
is used. The first step is to calculate the statistical proper- 
ties of the targets Rgg a nd P- The second step is to calculate 
the statistical property of the background clutter noise R^. 
The last steo is to calculate the filter coefficients by 
solving (3-3.15). They will be explained in more detail in the 
following . 



[Step I] Calculations of Statistical Properties of the Targets 
R ss and P: [14] 

The calculation of Rgg and P are relatively straight- 
forward. They depend on the type of target signal. Once the 
target signal of interest is given, as shown by several 
examples of point target and line targets in Fig. 3-3, 

Rgg and P are calculated by the following expressions: 
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3-3.3. Estimation 
Point at Lower 
Center Location 
of Search Box. 

Fig. 3-3#l,2 f 3. Three* Possible Estimation Point Locations in 
Search Box £>ut of 9 Possibilities in a 3x3 Search Box) . 























Respect to the 
Search Box in 
Fig. 3-3.1. 



Respect to the 
Search Box in 
Fig. 5-3.2. 



Respect to the 
Search Box in 
Fig. 3-3.5. 




Figure 3-3 Examples of Different Estimation Pixel Locations 
anu Different Target Types 
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(3-3-16) 
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where T = intensity of the point target. 

The second case is for a horizontal line target, as 
shown in Fig. 3-3. 7 (a), and for a central estimation pixel. 



Under these conditions, s 4 , s^ and s^ are not zero; all other 
s's are zero. 
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(3-3.20) 
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The amplitude of the line target s 4 , s ^ , s^ can be 
either deterministic or statistical. In many applications, a 
deterministic target amplitude is considered. If the interest 
is only in the location of the targets, s 4 = s^ = s^ = constant. 
If the amplitude of the target signal is known, s 4 , s 5 , and s 6 
will be proportional to its spatial variations. 
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Statistical amplitudes are useful in cases where a large 
number of targets are expected and their amplitudes are 
described statistically, for example, by a probability 
density function. This problem will be discussed further 
in Chapter V where a target detection procedure is devel- 
oped. 

Other types of targets can be handled in a similar 

manner . 



[Step II] Calculation of R^: 

Two procedures are used in this thesis to calculate 
the covariance matrix of the background clutter noise, 
depending on the type of images used in the analysis. 

1. Computer Generated Images: 

It is assumed that background clutters in real 
infrared images can be modeled by Markov processes. Their 
autocorrelation functions can be described by the following 
closed-form mathematical expressions: 

First order Markov process: 



Q Cm, n) = a c ^ p| n l 



( 3 - 5 . 21 ) 



Second order Markov process: 

Q(m,n) = cr^ p^ m l p^ n l cos (6 H m) cos (3 y n) (3-3.22) 

? 

where a~r. = variance of correlated noise which is used 
to model the clutter. 

p^ = horizontal correlation coefficient between 
two adjacent pixels. 
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P y = vertical correlation coefficient between two 
adjacent pixels. 

3 “m = ^ or i zonta l spatial frequency (for 2nd 
order Markov process) 

gy = oj n = vertical spatial frequency (for 2nd order 
Markov process) 



These two dimensional correlated random noise sources 
can be simulated by a recursive procedure described in Appendix 
A to generate simulated background clutter images. For these 
images, can be calculated easily [14], 

2. Real World Infrared Image: 

If real world infrared images are given, their sta- 
tistical characteristics must be calculated directly from the 
given images, using eq . (5-3.12) and then averaging it over a 

large number of pixels to obtain its expected value. It must 
be pointed out if the number of targets in the image is high, 
and their amplitudes are high, they will cause error in the 
calculation of the statistical properties of the clutter. 

The analysis of a real world infrared image was made by 
Captain D. Hilmers in his thesis research [50]. 



[Step III]. Design of Filter coefficients: 

After the calculations of R ss , P and R^, the filter 

coefficients are designed by solving eq. (3-3.15). 

5 . Figure of Merit for Clutter Suppression - 
Processing Gain 



To evaluate the performance of the designed filter, 
the processing gain (PG) was selected as a figure of merit. 
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The PG definition is: 



(Signal power/Noise power) At filter output 

PG = 

(Signal power/Noise power) At filter input 

( 3 - 3 . 23 ) 

The PG is given by the following equations: 

D 2 

Input signal power at a given pixel = s = Pg 

x ■ • • , D 2 2 1 

Input noise power at a given pixel = p 



N F N i 



'N 



is noise variance, defined by (3-3.22) 



is noise mean 



T,.„ 2 



Output signal power at estimation pixel = (S *W) " = Pg 
where s 1 was defined by eq . (2-2) and w by eq . (3-5.5). 

Output noise power at estimation pixel = 

e[(n t -w) 2 ] = e[w t -n-n^w] = w t -e[n-n t ]-w = w t -r nn -w = p N 

where R^ defined by (5-5.13). 



Therefore, the PG is: 



PG = 



(p s /p n ) 
0 0 



(s t -w) 2 /(w t .r nn -w) 



9 7? 

s“/(a N _ + y N “ ) 



(3-5.24) 



(Pg _ / p^j _ ) 

where, and were presented in Section II. E. 2 
For a point target and zero mean noise, the PG will become: 

l p • 
nN 



PG = 



s 2 • w 2 (MES , NES) / 0^ W T .R m p n .W w“(mES,NES) 



s 2 /(o n 2 ) 

p 



. R P . W 

k nn v 



(3-3.25) 



where R„„ = aj 1 • R^ and w(. MES, NES ) is the filter 
coefficient at the estimation pixel. 



4 . Results of Computer Simulation 

Properties of nonrecursive spatial MMSE filters 
and their performances in clutter suppression and target 
detection, when combined with the thresholding process, are 
investigated by computer simulation. Two types of images 
were used. One is a computer generated image containing 
spatially correlated noise represented by Markov models to 
simulate background clutter. The other is a real world 
infrared image. 

The search box size is limited to a 3x5 square. 

Point targets are used in most of the computer simulation. 

The estimation pixel is usually in the center of the 3x3 box. 
It is moved to the side of the 5x3 only when the estimation 
pixel is along the edges of the image. 

Investigation of the filters is carried out over a 
range of input conditions : 

SNRI = Input power signal to noise ratio 
= -21 db to 0 db 

= Correlation factor = 0.999 to 0.4. 

For target detection, a probability of false alarm (Pp A ) in 
the range of 10 ^ to 10 ^ is selected as a parameter. 

The computer simulation results are presented in 
Chapter VI where both single frame spatial filter results and 
multiple frame spatial - temporal filter results are presented 
together. Two types of presentations are used. The first 
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type is quantitative and does not use two dimensional pic- 
tures. Filter properties such as filter coefficients will 
be presented as functions of input signal to noise ratio and 
correlation factor. Filter performance, such as processing 
gain for clutter suppression and probability of detection 
for target detection, is related to input signal to noise 
ratio while the correlation factor and probability of false 
alarm are selected as parameters. The second type of presen- 
tation is pictorial and uses two dimensional images. Computer 
printouts using nine gray levels are used to present the 
unprocessed image, processed image after the application of 
spatial filter and images after further processing by thresh- 
olding . 

D. NONRECURSIVE SPATIAL MATCHED FILTERS 

1. Principle of Maximizing Signal to Noise Ratio [Ref.- 27-29] 
The objective of the spatial matched filter is to give 
a processed signal, Y = fy(m o ,n Q ) which will provide the max- 
imum output target signal to total noise ratio. The filter 
design procedure is presented in the following paragraphs and 
is an extension to two dimensions based on the theoretical 
derivation of Appelbaum [29J. 

[Step I] Unprocessed Image: 

The configuration of the nonrecursive spatial matched 
filter is the same as that of the MMSE filter. The description 
of unprocessed image and filters is the same also as presented 
in Section B . 1 . 



57 



(3-4.1) 



X = ( x 1 > x 2> . . -Xjj 



X = S + N 



(3-4.2) 



The filter is described by 



W = [w, w ? . . . w, ] 

1 y - j 



(3-4.3) 



[Step II] Filtered Image and Output Signal to Noise Ratio: 
The output of the filter is 



y = X *W 



(3-4.4) 



y contains two components: 



v = v + v 
S N 



where 



v 



2 = output signal voltage = S ♦ 



W 



T 



v^j = output noise voltage = N *W 



(3-4.5) 

(3-4.6) 

(3-4.7) 



The output signal power is | v ^ j “ . The output noise power is 

2 : 
n I 



the expected value of lv 



P n = E[ v n ’ v n ] = E [ | N T * W | 2 ] = E[W T *N*N T W] 

= W 1 E [N N ]W 

= W TA R nn .W (3-4.8) 

(the superscript * refers to the complex conjugate wherever a 
superscripted variable is complex; otherwise the * can be ignored) 



where 



R..,, is the covariance matrix of the noise given bv 
NN ° 



R 



NN 



E(N 1 N ] _) E(N x N 2 ) . . . E(N 1 N I ) 

ecn 2 n 1 ) e(n 2 n 2 ) . . . e(n 7 n I} 

E(N i N 1 ) E (NjNj) 



(3-4.9) 
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assuming N has a zero mean. The output power signal to noise 
ratio is: 



SNR 



out 




w"^*R *W 
k NN 



(3-4.10) 



The matched filter coefficients can be designed by maximizing 
(3-4.10). However, a more effective way to derive the matched 
filter design equation is to diagonalize the matrix R,^. 



[Step III] Diagonalize Output Signal to Noise Ratio: 

Let A = the linear unitary transformation which 
diagonalizes the matrix R^ which is equivalent to decorre- 
lating the input clutter noise, i.e., 



A*R nn A T = I and R nn X = (A T -A*) (3-4.11) 

The new filter configuration is shown in Fig. 3-4. .The 
resultant transformed input signal and noise are: 



S' = A-S 
N’ = A* N 



(3-4.12) 

(3-4.13) 



Let W' be the new transformed filter, then the output signal 
and noise voltages are 



v s = W' T .S' = w ,t *a*s 


(5-4.14) 


V N = W ,T -N’ = w ,t -a*n 


(3-4.15) 


Comparison of (3-4.13) and (3-4.14) with (3-4.6) 
the relation between W and W' can be found as 


and (3-4.7) , 


W = A T -W = W' a t 


(3-4.16) 
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Figure 3-4 Original and Transformed Filter Con- 
figurations Used in the Derivation of 
Spatial Matched Filter Design Equation 
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The output noise power after the transformation is 



P n = E[|W' T #N'| 2 ] = E[W’* T N'* N ,T W'] 

= W' * T E[N’ * N ,T ] W’ (3-4.17) 

Since the transformation A diagonalized the R^ and decorrelates 
the noise 



E[N'* T N'] = I 



(3-4.18) 



and 



p = W’ * T W 1 
n 



(3-4.19) 



The output signal power after the transformation is 

, . , 7 



v r 



(3-4.20) 



Therefore, the output signal to noise ratio is 

2 



v r 



SNR 



out 



W' * T W 



(3-4.21) 



[Step IV] Equation of Matched Filter Design: 

The design equation for the spatial matched filter 
is obtained by finding the maximum value for (SNR) out . To 
accomplish this, the Cauchy Schwarz inequality is applied to 
[5-4.14) . 



where 



v. 



< |W 



T , 2 , ,2 

i 1 r Cl ^ 



< I 2 = w 1 * T W’ 



I I 2 = | n , *T 



S' 

, 2 



(SNR) 



ou t ^ i i 



< IS 



, . 2 



(3-4. 22) 
(3-4.23) 
(3-4.24) 
(5-4.25) 
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Therefore, the maximum output SNR is 



CSNR) = IS 
^ 'max • 



2 



(3-4. 26) 

It will be shown in the following that the special 



filter 



W = cS ' * 



(3-4.27) 



where c = an arbitrarily selected constant 



will provide the maximum (SNR) 0U - t p Ul - • For this filter 
v g = cS ' T S ' = c| S ’ | 2 



(SNR) 



Pn 



out 



= c 2 |S' | 2 



v r 



= I S' I 2 = (SNR) 



n 



max . 



Thus, it has been shown that the special filter, W' = cS 1 , 
is the matched filter after the transformation of A. The 
matched filter itself is expressed by: 

W. fC = A T W = A T • c • S ' * = A T -c A*S* 

Mr 

= c A T -A*S a 

Using the definition of (3-4.11) 




(3-4.28) 

This is the design equation for the nonrecursive spatial matched 
filter. Using this matched filter, the maximum signal to 
noise ratio is: 



SNR = S T R ' 1 S* 
max NN 



(3-4.29) 
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It is worthwhile to note that the maximum output signal to 
noise ratio depends only on the target signal and clutter 
noise characteristics. In active systems, such as radar and 
communication systems, the signals can be designed to yield 
very high SNR. On the other hand, in passive systems, such 
as passive infrared and optical surveillance systems, the 
target signal can not be designed. The maximum SNR is given 
by (3-4.29). 

2 . Design Procedure of Spatial Matched Filter 

To design a nonrecursive spatial matched filter, a 
procedure somewhat similar to that of the nonrecursive spatial 
MMSE filter is used. This is readily clear by comparing their 
design equations. 

For MMSE filter: W wp = (R + R ) _1 P (3-3.15) 

(Wiener filter) 

'For matched filter: W,. p = cR^ ^•S* (3-4.23) 

A three step procedure is also used. 

[Step I] Selection of c and Calculation of S* : 

Although the selection of c is arbitrary from the 
viewpoint of matched filter design, its selection is important 
since its value will affect the threshold setting for target 
detection. This will be discussed in more detail in Chapter VI. 
In general, the value of c is chosen so that the filter 
coefficient corresponding to the estimation pixel is normal- 
ized to unity. 
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The determination of S* is 
will be given. For a point target 

S* T =[0000S000 

For a horizontal line target as in 



very simple. Two 
as in Fig . 3-3.1. 

OJ. 

Fig. 3-3.7 (a) . 



examples 



S* T =[000SSS000] 
where S = amplitude of target 



Other types of targets can be handled in similar ways. 



[Step II] Calculation of : 

This is exactly the same as described in section B.2. 

[Step III] Calculation of Filter Coefficients: 

The filter coefficients are calculated by solving 
equation (3-4.28) . 

3 . Processing Gain 

The processing gain (PG) of the matched filter is 
derived in a way similar to the derivation in section B.3., 
using equation (3-3.24) as a generalized PG equation. For 
a point target, equation (3-5.25) can be used to determine 
the PG. When the constant c in the matched filter is selected 
as suggested in D.2., the PG of point target will simply 
become : 



PG 




1 




W 
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For a point target, assuming noise with zero mean, the P.G. 
can be shown to be independent of input SNR. Defining 
S = s[0 0 . . 1 . . 0] = s*S j. and = — 2 an< ^ 

substituting these definitions into eq. (3-4.29) yields: 



PG 



SNR 

SNR 



(output) 
( input) 



SNR 

SNR 



max 

( input) 



2 c T D 

s S I * R NN 

2 , 2 
s /a 



-1 



S 



T 

I 




S 



I 



This last expression of P.G. is independent of signal and 
noise power and depending only on the correlation factor of 
the correlated noise and the target location in the search 
box . 

4 . Results of Computer Simulation 

Properties of nonrecursive spatial matched filters 
and their performance in clutter suppression and target 
detection when combined with the thresholding process were 
also investigated using computer simulation. 

The items of investigation in terms of search box 
size, background clutter characteristics, target types, filter 
coefficients, performance figures of merits, are the same as 
in cases of nonrecursive spatial MMSE filters. They have been 
described in Section C.4. and will not be repeated here. 



E . SUMMARY 

Two types of statistical spatial filters were discussed 
in this chapter: the spatial MMSE filter and the spatial 
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matched filter. Both filters offer some signal processing 
advantages . 

The MMSE filter, designed to minimize the mean square 
error, will offer the best spatial estimate of a given target. 

The matched filter, designed to maximize output signal to 
noise ratio, will result in the best enhancement of SNR suit- 
able for detection purposes. 

It turns out that for the detection of a point target, 
both the MMSE and the matched filter offer the same processing 
gain. This will be shown in Chapter VI. 

In general, the type of filter to be selected will be 
determined by the practical processing requirements for a 
given problem. For target enhancement and restoration tech- 
niques which require accurate estimation of the target's shape 
and intensity, the Wiener filter is preferred. For target 
detection requirements which allow distortion of target's 
shape and intensity, the matched filter is preferred. 
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IV. MULTIPLE FRAMES FOCAL PLANE PROCESSING FOR 

CLUTTER SUPPRESSION - NONRECIJRSIVE STATISTICAL 
TEMPORAL FILTER AND SPATIAL - TEMPORAL FILTERS 



A. INTRODUCTION 



If an imaging system can be stabilized or tracked to 
follow a scene, several successive frames of images can be 
used for various image processing purposes. Signals from 
multiple frames of images can be used in several ways. The 
simplest way is to use one pixel from each frame; signals 
of the same pixel from successive frames are used together 
in what amounts to temporal processing. A more involved 
way is to use several pixels from each frame; for example, 
signals from nine pixels in a 3x3 search box as developed 
in Chapter III are used together with signals of the same 
search box location from consecutive frames in what amounts 
to a combined spatial - temporal processing. 

In this thesis, an intermediate process is developed. A 
spatial filter is first used to process several successive 
frames of images individually. Then, single pixels from 
each of these frames already processed by the spatial filters 
are used together to perform a temporal filtering. The 
spatial filters are designed by statistical criteria as pre- 
sented in Chapter III. The temporal filter will be designed 
in a similar manner, using the same principle - that if the 
temporal noise is correlated, it can be suppressed by a 



nonrecursive temporal filter designed either by a minimiza- 
tion of mean square error criterion or by a maximization of 
signal to noise ratio criterion. The objective of these 
temporal filters is the same as other focal plane processing 
techniques developed in this thesis, namely, signals from 
single frame and/or multiple frames of images are processed 
to suppress the background clutter so that weak targets buried 
in them can be pulled above the noise, thus enabling a thresh- 
old to be applied to initiate the target detection process. 

In Section B, a class of nonstat is tical , multiple frame 
focal plane processing techniques, called "frame differ- 
encing" techniques, are first presented. An analysis is then 
made of their effects on the statistical properties of the 
images. These nons tatist ical techniques will be compared with 
the performance of statistical temporal filters developed in 
this thesis. In Section C, three dimensional statistical 
filters using multiple frames of images are developed. Spe- 
cifically, a sequential combination of several spatial filters 
and one temporal filter is developed. In Section D, an accum- 
ulative implementation of the nonrecursive temporal filter is 
suggested which may simplify the memory requirements in hard- 
ware implementation. In Section E, computer procedures for 
investigation of these multiple frame, focal plane processing 
algorithms are described for both simulated images with pre- 
scribed statistical properties and real world images without 
a priori knowledge of their statistical properties. 
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However, the simulated results will not be presented in 
this chapter. They will be presented in Chapter VI, which 
includes the simulated results for both the single frame 
spatial filters and the multiple frame spatial- temporal 
filters . 

B. NONSTATISTICAL TEMPORAL FILTERS - FRAME DIFFERENCING 
1 . Introduction 

Multiple frame processing has been actively developed 
in the past few years for space surveillance applications, 
especially for those systems in which the image sensors 
platform is stabilized to cover the same area for certain 
periods of time. These systems are generally known as the 
staring, or step- staring , systems. In these systems, there 
is a great deal of correlation of background clutter among 
consecutive frames of images. A simple subtraction of two 
neighboring frames should eliminate completely, in principle, 
all background clutter. This is most effective in revealing 
those weak targets which are moving and do not dwell in the 
same pixel between two frames. However, in real world situ- 
ations, sensor systems cannot be perfectly stabilized. Also, 
not all background clutter is stationary. An example is 
clouds which, move. Consequently, simple subtraction of 
two neighboring frames is not adequate in many practical 
applications. More involved subtraction processes using 
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several frames of image have been developed. They are 
commonly known as the "Frame Differencing" methods [10]. 



"Single Differencing" using two frames = 



F 



1 



F 



2 



"Double Differencing" using three frames: 



F i * 2F, + F- ='(F 1 - F 2 ) - (F 2 - F 3 ) 



"Third Differencing" using four frames: 




"Fourth Differencing" using five frames: 



etc . 




where F^ refers to the signal from the same pixel from the 
"i"th frame. 



from the simple concept that "non-moving" clutter can be 
cancelled out by frame to frame subtraction. The weighting 
factors for individual frames are fixed and are independent 
of the specific nature of the images. Since the statistical 
temporal filters to be developed in this chapter are designed 
according to the specific statistical properties of a given 
set of images, a statistical analysis of the standard "frame 
differencing” techniques will be made so that they can be 
compared with the statistical temporal filters developed in 
this chapter. 



These "frame differencing" methods were developed 
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2 . Statistical Analysis of First Differencing Technique 
Since the same pixel location is used in "frame 
differencing" techniques, its coordinate (m,n) will not be 
used in the following analysis. Instead, a number k is used 
to describe the "k"th frame. In the following analysis, 
special attention is given to the jitter problem which is 
modeled by a random noise component from frame to frame, i.e., 



a, 8 are two constants 

nC0,l) = a Gaussian noise component of zero mean 
and unity variance 

x(k) = background clutter level at a given 
pixel in the k-th frame 

x(k+l) = background clutter level at that pixel 
in the (k+l)-th frame 



In eq . (4-1) it is assumed that the noise intensity at a pixel 
in the (k+1) frame is proportional to the intensity in the 
same pixel at frame (k) with proportionally constant a. The 
second term in eq. (4-1) represents the noise added to that 
pixel due to the sensor jitter. 



The constants a, 8 can be rewritten in different forms 



compatible with the notation used in the Markov model in 
Appendix A. 



x(k+l) = ax(k) + 3n(0,l) 



(4-1) 



where 



x(k+l) = p K xCk) + a K /l- p K 2 n(0,!) 



(4-2) 



where 



a = = temporal correlation factor 



7 

3 = a K /l- p K “ 

Oj- = temporal standard deviation. 



71 



In the following, the statistical properties of the 



first differencing method will be derived. 

Let z Ck) = first differenced signal 

= x (k) - x (k+1) 

= (1-P K ) x (k) - a K /l-p K 2 * n(Ofl) (4-3) 

The mean of z(k) is zero since x(k) = x(k-l), (4-4) 

where the symbol x is the expected value of x,E(x), taken 
by averaging over the two dimensional image. The variance 
of z(k) is 

E [ z 2 (k) ] = E [ ( 1- p K ) 2 x(k) 2 + ^ k 2 (1-P k 2 ) n 2 (0,1) 

+ 2 (l-p K )x (k)a K /l-p K 2 n (0 , 1) ] 

i l i 2 2 

= (l'Pj() -a x~ + a K' , *' 1 ~ p K*"^ a n~ + E[X(K)n(0,l)] ‘constant 



The last term is zero because x(k) and n(0,l) are assumed to 
be uncorrelated and since n(0,l) has zero mean. 

E[z”(k)] - a x (1-pjr) + o ” ( 1 - pjr ” ) p n 



Since the variance of n(0,l) is unity, 

E [t 2 (k) J = a x 2 (l-p K ) 2 + a K 2 ^ 1 ' p K 2 ^ ) 



where 2 . 

= temporal variance 

o y ~ = variance in one frame 



when 



Q O 

E[t 2 OOJ = a 2 wCl-p) 



(4-5) 



(4-6) 
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3 . Figure of Merit for Frame Differencing - 
Processing Gain 

A figure of merit, "processing gain," is defined to 
describe the performance of frame differencing. The same 
definition is used for processing gain for all filter tech- 
niques developed in this thesis. 



where 



PG = 



(S/N) 



0 



(S/N) . 



2 



(S/N) 



2 D 



0 



= output 



(S/N)/' B input p 



power 

ower 



s ignal 
signal 



to noise ratio 
to noise ratio. 



(4-7) 



The "signals" considered in this thesis are point target signals. 

In this chapter, it is assumed that these point target signals 

move at least one pixel from one frame to the next. 

2 2 
Let = signal power in the first frame pixel = S“ 

S 7 “ = signal power in the same pixel in the second 
frame 

= 0 because the target has moved away 

. ? 2^7 

•• S“ ut = output signal power = - S ? " = S“ 



, 2 

'input 



S7 ^ = input signal power = S^“ = S' 



The input noise power is 



N/ 



- a 



X + ^X 



where = input noise variance 

* 2 

= input noise mean power 
The output noise power is 

N n 2 = E[z 2 ] = a Y 2 (l-p Y ) 2 + o K 2 (l-p Y ) 2 



'0 J U X ^ x ” K X 

with zero mean as shown by (4-5) 



K ; 
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Therefore, the processing gain is 

s 2 



PG = 



(S/N) 2 q _ a x 2 Ci- p K ) 2 ^ a^(l-p 2 ) 



(S/N) 2 . 



2 2 
CT X + U X 



PG 



i*Y 2 

1 + 

°x 



? ° Y "* ? 

<>^' + ^ ^-°K > 



(4-8) 



For the special case of 



'K X 



? 



PG = 



1 + C<r) 



2(1-P k ) 



(4-9) 



This "first differencing" processing gain will be used in 
section C.2 when it is compared with the processing gain of 
statistical spatial - temporal filters. 



C. STATISTICAL SEQUENTIAL SPATIAL -TEMPORAL FILTERS 

In this section, a three dimensional statistical filter 
is developed. It consists of a sequential application of 
several spatial filters to successive frames of images fol- 
lowed by a temporal filter using these processed images. In 
contrast, another three dimensional filter processing signals 
from several frames of image all at once will be presented in 
the next section, D. 
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1 . Filter Configuration 

The configuration of the sequential spatial- temporal 
filter is presented in Fig. 4-1. Its operation is described 
as follows: K frames of images consist of point targets 

contaminated by background noise. 

Step I - Spatial Filters for K Frames : 

The detailed design procedure for spatial filters have 
been presented in Chapter III. The highlights of their 
operation will be repeated here for completeness of the 
description of sequential spatial- temporal filters. 

[Unprocessed Signal] 

A two dimensional array of unprocessed signals in a 
selected search box around the estimation pixel (m,n) is 
taken from each frame. 



This two dimensional array of data is represented by a 
one dimensional vector. 



X(m,n) = Xj 




(4-10) 



*— X NSB , 1 , X NSB,MSB J 




. . X . . . 

1 




(4-11) 



j = MMAX(n-l) + m 



n = 1 , . . . NMAX 



m = 1 , . . . MMAX 
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Three Frame Sequential Spatial-Temporal Filter 



Spatial 
Filter W 





Temporal 
Filter W 





Spatial-Temporal 
Filtered Output 

▼ z 



Figure 4-1 Configuration of Nonadaptive Three Dimensional 
Nonrecursive Sequential Spatial -Temporal Filter 
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i = MSB • (An- 1) + Am 
An = 1 , . . . ,NSB 
Am = 1 , . . . , MS B 
I = NSB-MSB 

as presented in Fig. 4-2. 

The unprocessed signal X contains target signals 

S = [s 1> s 2> • • • >Sj] (4-12) 

and clutter noise 

N = [n-j. n 2> • • • >n j ] (4-13) 

[Spatial Filter] 

A linear two dimensional filter is described by a two 
dimensional array of weighting coefficients 



W (m , n) 



w ll, w 12, ’ ‘ ' ,W 1,MSB 
w 2 1 , w 2 2 , • • • ’ w 2 ,MSB 



_ W NSB,1 



, w 



MSB, MSB. 



(4-14) 



It can be described also by a one dimensional vector 



W = W(i) = [w x w 2 ... ,w T ] (4-15) 

These filter coefficients can be designed either by the 
minimization of mean square error criteria or by the maximi- 
zation of signal to noise ratio criteria. 



[Spatial- Filter Output] 

The output of the spatial filter is 
y (m , n) = y = W* X^ 



(4-16) 
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The Image Configuration 
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T 


1 
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X (m , n ) = S (m , n ) 

+ N(m,n) 



Image Size 
(MMAX • NMAX) 




Figure 4-2(a) Coordinate System for the Image 

(b) Coordinate System for the Spatial Filter 
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Step II - Temporal Filter: 



[Input Signal] 

Filtered signals at the same pixel location (m,n) of K 
frames are put together to form the input signal to the 
temporal filter. 



Y 





(4-17) 



[Temporal Filter] 

A linear one dimensional filter is described by a one 
dimensional vector. 



W K = W K;( m > n ) = f w Kl , W K2 , ’ ' ' ,w KK^ (4-18) 

The filter coefficients are designed by the same procedures 
used for spatial filters by following either the minimization 
of mean square error criterion or the maximization of signal 
to noise ratio criterion. 



[Combined Spatial-Temporal Filter Output] 

The final processed signal at pixel (m,n) after the 
sequential application of K spatial filters and one temporal 
filter using K frames is described by 

z(m,n) = Y*W k T (4-19) 

2 . Statistical Analysis of Two Frame Sequential 

Spatial -Temporal Filter 

A statistical analysis of the two frame sequential 
spatial- temporal filter will be analyzed first so that it 
can be compared with the result of the nonstatistical first 
differencing method analyzed in Section B. 



79 



The analysis o£ the more general K frame sequential spatial- 

temporal filters will be presented in the next section. 

a. Statistical Model of Two Frame 
Single Differencing: 

The outputs from the spatial filtering are 
described by: 



y (1) = X(1)*W^ where X(l) denotes X of frame 

y(2) = X(2)*W^ where X(2) denotes X of frame 

In matrix form, they are 















y (D 




w 


0 




X T (1) 


y (2) 


— 


0 


w 




X T ( 2) 















CD 

( 2 ) 



(4-20) 



To calculate the noise power of the final output, the 
covariance of the spatial filtered image, Ryy> is needed. 
Using (4-20) : 

Ryy = E [Y T * Y] 



* 
















— 


w 


0 




X T (1) 




X(l) 




T 

W 1 


0 


0 


w 




X T ( 2) 




X (2) 

— 




0 


W T 





w 


0 




R XX (DD 


R XX (D2) 




"w T 


0 




0 


w 




R XX C2,D 


r xx (2,2) 




0 


W T 



where R.yx (i , j ) = E [X^ CDX (j ) ] = covariance of [X] 



(4-21) 

(4-22) 



It is assumed that both X and Y have zero mean and that frame 
(1) and frame (2) statistical properties are the same. 
Therefore : 



R xx (l,l) = %X ( 2 ’ 2 ) = R 



(4-23) 
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I 



Using the first order Markov model, as in section 
IV. B. 2, ((4-2) the temporal statistical characteristics are: 



X T C2) = p k I-X T ( 1) + I-a K /T^T K 2 n(0,l) (4-24) 

where I is the unity matrix. 

Multiplying (4-24) by X(l) and taking the expected value 
yields : 

E [X T ( 2) • X ( 1) ] = E { [p K I X T (1) + I a^T*^ 2 n (0 , 1) ] X( 1) } 



R 



xx (2,l) = P K *E [X (l)X(l) ] + I V^Pk" n(0,l)X(l) 

(4-25) 

The second term is zero since n is Gaussian with zero 
mean and is not correlated with X(l). Therefore, 

R xx (2,l) = R xx (1,2) = P k E[X T (1)X(1)] = (4-26) 

Substituting back into (4-21), 



R 



YY 



w 


0 




R xx 


p k r xx 




r t 

w 1 


* 

0 


0 


w 




1 TD 

pa 


R xx. 




0 


w - 



"w R xx W T 


P K W R XX wT " 


_ P K W R XX W 


w r xx w t 



(4-27) 



T ? 

W R xx W = Oy is the noise variance of spatial 



filtered images since 
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Y 

Therefore , 
R 



E [y • yj = E[ CWX T ) CWX T )] = E[W-X T *X-W T ] 

w-e[x t -x]*w t = w R x y 



YY 





- - 

a Y 


p K a Y~' 


2 


■1 


= 


2 


2 


- CTy » 






- p K a Y 


a Y J 




_ p k y 



(4-28) 



b. Statistical Analysis of Two Frame 
Single Differencing 

The temporal filter used in this analysis is the 
simple first differencing operation, so that the temporal 
filtered output is: 

w K = [1,-1] 
z = [1,-1] 



yd) 

y(2) 



= y(i) - y(2) = w K y 



The variance of z is: 



z 2 = E [ z 2 ] = E[W k .Y T * Y*W k T ] = W k E[Y T *Y]W k T 



(4-29) 



(4-30) 



Substituting (4-28) into (4-30) will yield: 

1 P K 



°z = W K‘ R YY' W K T = [d-l]a. 2 



Y 



PK 



1 

-1 



O Y 2(1- Pjr) 



(4-51) 



It can be shown that the same result is true for 
a sequential filter by applying first the temporal differencing 
followed by the spatial filtering. 
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c. Processing Gain 

It will be assumed again that the target signals are 
point targets which move at least one pixel from one frame 
to the next, so that 

s Cl) = s 
sC2) =0 
as in IV. B. 3. 

Therefore , 

s 



b . 

in 



out 



s ( 1 ) - s ( 2) = s 



The input noise power is given by 
2 2 

N. = a v assuming X has zero mean, 
in a & 

The output noise power is given by (4-31). 

Finally, the processing gain is: 

s 2 






PG =■ 



a y “ 2 (l-p 



K 



X 



a Y “2(l-p K ) 

a x' 



'X 

f Y 



a y 2 2(1-p k ) 



W R JL W T 



XX 






= (PG of spatial filter) • (PG of first 
differencing) 



temporal 

(4-32) 



using the relation R 



‘XX 



= a 



X 



R P 
R XX‘ 



(4-32) shows that the processing gain of the sequential 
spatial- temporal differencing is a product of two individual 
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processing gains whereas the temporal differencing alone 
possesses only the PG related to the second term in eq. (4-32). 

3 . Statistical Analysis of K Frame Sequential 

Spatial-Temporal Filter 

In this section, a more general K frame sequential 
spatial- temporal filter is analysed. In this case, the 
temporal filter used is designed by a statistical process 
instead of the nonstatistical "first differencing" temporal 
filter used in the previous section. 

a. Statistical Analysis of Spatial Filters 
of Different Frames 

The input temporal signal from K successive frames 
is described by spatially filtered output Y where: 



Y = [y I 7 7 • • 



Yk ] 



which results from applying spatial filters to K frames of 
unprocessed signal X. 



y(i) 

y(2) 



y(K) 



W 0 0 0 X T (1) 

0 W O' 




w 

W X T (K) 



(4-53) 



or, 



T ~ ~T 
Y 1 = W X 1 



C4-34) 



Again, the covariance of Y is needed for calculating the noise 
power of the images processed by the spatial filters. 
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If we assume Y has zero mean, 



R YY 

Covariance [Y] = E[Y T -Y] = E[W-X T -X-W T ] 

= W*E[X T *X]W T (4-35) 

■ s xx ftT 



The same first order Markov model is used for the correlation 
of two successive frames. 

X T (K) = p K I X T (K-1) + I ^ K ^U~ r P K 2 n(0,l) 

Therefore, similar to eqs . (4-Z5) and (4-26), 

E[X T (K) -X(K-l)] = p K R xx (4-36) 

Similarly, it can be shown that 

E[X T (K) X (K- 2) ] = p K 2 R xx 



In general. 



E[X T (K) X (K-j)] = P K J R xx 



(4-37) 



Substituting (4-37) into (4-55) and using the following 
two relations - 



R YY = W • 



R 



XX 


p k r xx p k R XX • • 


p K-l 
K 


R xx 


D 

K XX 


R xx p k r xx 


pK-2 

K 


R xx 


f 1 R_ 









and 



XX 



W 



T 



(4-38) 



= W’R XX W = variance of spatial filtered image. 
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1 



R yy ~ a y 



K 



'K 



iK* 1 



or 

R 



YY 



2 

^ Y ^YY 



'K • 



,K-1 

K 

K-2 



( 4 - 38 ) 



b. Statistical Analysis of Temporal Filters 

Similar to the spatial filter cases in Chapter 
III, two statistical criteria can be applied to design the 
temporal filter. 

If the minimization of mean square error criterion 

is used, 



W K = [ R yy„] P 



YYK J 



K 



( 4 - 39 ) 



where R yyK - R NNK + R SSK 

R NNK = covar i ance [^}(] 

R$sk = covariance [S^.] 

P K = E[Y-S k ] 

As before, target signal and clutter noise are 
assumed to be uncorrelated. 

If the maximization of signal to noise ratio 
criterion is used, 

_1 - ( 4 - 40 ) 



W K C ^ R NNK^ S K* 
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where C = selected constant 

= target signal vector in the same pixel along 



K frames 

where k = 1 , 2 , . . . K 

The final output of the sequential spatial - temporal filter 



[S K1, S K2, ’ * s Kk, • ,S KK ] 



is 



z = Y* W 



K 



The noise output power of z is given by 

z = w k ,r yy* w k = w k ,r yy* w k a Y 



(4-41) 



c. Processing Gain 

Assuming again that the target signals move at 
least one pixel from one frame to the next - 



s . = W„»S„ 
out K K 



where = [ 0 , 0 , . . . s , 000 . . ] for a moving point target 
In most applications of point target detection, 



out 
2 



S = s. 

in 



N. 

in 



N 



out 



= one frame input noise power = o 

A 

= output noise power after multiple frame 
filtering = a ^ 



PG = 



CS/N)' 



out 



X 



(S/N) 2 



(4-42) 



m 



'X 
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where o. 



Y 



W K* R YY 



• W 



T 

K 



a 



2 

Y 



• W 



K* 



R 



Pk 

YY 



W 



T 

K 



W. R xx W T = a x 2 -W-R xx -W T 

CT X t w * R xx' W ^ w k* r yy‘ w k ^ 



Finally, 



D. 



ru - y * T- t 1 . 4-4 

CW.Rx^.W 1 ) ( w k‘ r yy' W K } 

= (PG of spatial filter) (PG of temporal filter) 

AN EFFECTIVE WAY OF HARDWARE IMPLEMENTATION 
OF SEQUENTIAL SPATIAL- TEMPORAL FILTER 



In the previous section, a special class of a three 
dimensional spatial - temporal filter was developed in which 
nonrecursive spatial filtering and nonrecurs ive temporal filter 
were carried out sequentially. They are described by the 
following operations. 



y CD 

y (2) 

yCK) 



i 

= E x(i,l)w(i,l) 
i= 1 

= E x(i,2)w(i,2) 
i=l 



= .| i xCi,K)w(i,K) 
K 

y(k) (k) 



= X(1) 
= X(2) 



= X (K) 
T 

Y W/ 



W 



W J 



w 



Spatial 



Temporal 



In general form, 

K MSB NSB 

z = Z w^(k) Z E w(m,n)x(m,n) 
v-i K m= 1 n=l 



S8 



This type of nonrecursive processing requires the memory of 
all K frames of image data. 

It might be possible to implement this type of three 
dimensional filter by only using one frame of memory and 
updating its signals in a recursive manner as follows: 

MSB MSB 

z (k) = z(k-l) + w„(k) Z Z w(m,n)x(m,n) 

* m=l n= 1 

This approach of combining a recursive implementation of 
nonrecursive temporal filter with nonrecursive spatial filter 
could be used when the temporal filter is designed by other 
criteria, such as by the Kalman filter approach, instead of 
the MMSE filter and matched filter approaches developed in 
this thesis. This recursive approach is suggested as a method 
to be developed in possible future work. 

E. COMPUTER SIMULATION PROCEDURE FOR EVALUATING 

SEQUENTIAL SPATIAL-TEMPORAL FILTERS 

The three dimensional sequential spatial - temporal filters 
have been evaluated by computer simulation using two types of 
images : 

a. Simulated images using Markov model clutter noise 

b. A real world infrared image. 

The results will be presented in detail in Chapter VI. In 
the following, software procedures for evaluation of these 
two types of images will be presented. 
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1 . Markov Model Simulated Images 



The following steps are used. 

a. Generate a two dimensional stochastic field with 
specified statistical properties to simulate the background 
clutter [Appendix A) . 

b. Add a set of targets in the clutter to form an 

image . 

c. Generate K frames of images with specified temporal 
statistical properties. 

d. Design spatial filter and then temporal 
filters according to either MMSE or matched filter criteria, 
optimized to the assumed clutter statistics and target char- 
acteristics . 

e. Apply sequentially the spatial filter and temporal 
filter to suppress the clutter and to enhance the targets and 
evaluate the processing gain. 

f. Calculate a threshold setting according to a spec- 
ified probability of false alarm assuming that the probability 
density function of the filtered image is Gaussian. 

g. Apply the threshold to separate targets from 
residual clutter. 

2 . Real World Infrared Images of Background Clutter 

Because the statistical properties are inherent in 
the images, no computer simulated clutter is needed. The 
following steps are used. 
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a. Read the k-th frame, calculate its statistical 
spatial properties such as mean, variance, probability density 
function and covariance matrix. The exact calculation method 
is presented in Ref. [30] . 

b. Design statistical spatial filters according to 
either the MMSE filter criterion or the matched filter criterion 

c. Add a set of targets to the given infrared images. 

d. Apply the spatial filter, evaluate the processing 
gain, and store the processed k-th frame. The detailed pro- 
cedure can be found in Ref. [30]. 

e. Repeat steps (a) to (d) for all the K frames. 

f. Calculate the statistical temporal properties of 
selected pixel locations from all K frames of spatially 
filtered data. 

g. Design a statistical temporal filter according to 
either the MMSE filter criterion or the matched filter criterion 

h. Apply the temporal filter, calculate the mean and 
the variance of the filtered image and evaluate the processing 
gain. 

i. Calculate a threshold setting for a specified 
probability of false alarm using the method to be described 
in Chapter V because the probability density function is 
generally non-Gaussian. 

j. Apply the threshold to separate targets, false 
alarms from residual clutter noise. 
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F . SUMMARY 



Three dimensional focal plane processing methods are very 
effective in background clutter suppression because the pro- 
cessing gains of spatial and temporal filters are multiplica- 
tive. A combination of spatial and temporal filtering should 
always improve the total processing gain. Furthermore, if 
the temporal clutter noise is highly correlated, which usually 
is true if the drift and/or jitter between successive frames 
is small, the processing gain of a temporal filter is usually 
higher than that of the spatial filters. An example will be 
presented in Section F of Chapter VI. 

In this chapter, the configuration, design procedure and 
performance analysis of nonstatistical "frame differencing" 
temporal filters developed elsewhere [10] were first presented, 
followed by the sequential statistical spatial- temporal filters 
developed in this thesis. It was found that the sequential 
spatial temporal filtering technique, when applied on real IR 
images, offers indeed superior processing gain compared to 
either spatial filtering or frame differencing. 
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V. TARGET DETECTION IN FOCAL PLANE PROCESSING 



A. INTRODUCTION 

In most infrared images, targets are deeply buried in 
noise which consists of strong background clutter noise and 
other random noise. The targets can be as weak as 70 to 80 
db below clutter. Predetection processing must be used to 
suppress the clutter to the point that "target" signal to 
"clutter" noise ratio is adequate for target detection. This 
was discussed in Chapters III and IV. 

Target detection usually consists of several steps. The 
first step is "thresholding" which is the subject of this 
chapter. This is used to separate targets from noise remain- 
ing after the predetection processing. Since both target 
signals and noise after filtering are statistical quantities, 
it is generally difficult to clearly separate targets from 
noise except in rare occasions when the signal to noise ratio 
is very high. In most cases, the intensities of targets and 
noise overlap. Any selection of threshold will miss some 
targets and include some noise (false alarms) . There is 
always a compromise between a complete detection of targets 
by using a low threshold at the expense of including many 
false alarms, and a complete elimination of false alarms by 
using a high threshold but at the expense of missing many of 
the targets. Two figures of merit are normally defined to 
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describe the false alarm and target detection. They are the 
probability of false alarm CPfa) an< ^ probability of detection 

(P D ) • 

A target detection process involves five types of related 
information : 

1. The PDF of filtered targets and noise (PDF is probability 
density function) 

2. Signal to noise ratio - SNR 

3. Threshold level 




Their relationships to detection have been very well developed 
in radar, sonar and elsewhere [31-33]. In this chapter, the 
one dimensional detection procedure is extended in a straight- 
forward way to two dimensions for application to focal plane 
processing problems. 

B. PREDETECTION FOCAL PLANE PROCESSING FOR CLUTTER SUPPRESSION 
Clutter suppression is one of the most important processing 
steps before the thresholding for target detection. Statis- 
tical spatial filter and spatial - temporal filters have been 
developed in Chapters III and IV respectively. Their impor- 
tant steps and statistical descriptions will be reviewed and 
summarized in this section. 

1 . Statistical Properties of Unprocessed Image 

The unprocessed image data taken from a two dimensional 
detector array of Am x An pixels are described by: 
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X 1 1 , X 2 1 ’ * • X MSB,1 

X = 

X 1,NSB * ' * X MSB , NSB 

In this study, this two dimensional array is written as a 
one dimensional vector for ease of programming in processing 
software . 

X = X (i) = [x lj x 2j . . . Xj] (5-1) 

where i = MSB’ (An- 1) + Am Am=l , 2 , . . . .MSB 

An=l, 2, NSB 

This transformation was previously presented in Fig. 3-2. 

The unprocessed image is considered to be the sum of target 
signals 

S = S(i) = [s, s 7 . . . , s T ] (5-2) 

and background clutter noise 

N=N(i)=[n,n ? ...,n T ] (5-3) 

Target signals can be either deterministic or 
statistical. It is assumed that their spatial shapes are 
deterministic such as point target orline targets of various 
orientations. Their intensity amplitudes may be deterministic 
or statistical with Gaussian probability density functions 
assumed. Deterministic targets are described by the vector - 

S = t s l, s 2 • • ’ S I^ 
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Statistical targets are described by several 



statistical characteristics: 



Mean vector = y<, = 


^ y Sl, y S2, * ' * y SI^ 


(5-4) 


Covariance matrix = 


R ss = E [ (S T - u s T ) (S - y s )j 


(5-5) 


Variance = Og“ 




(5-6) 


Autocorrelation function = Qg 


(5-7) 



The background clutter is statistical and described 



by the following characteristics: 

Mean vector = = [u N1 u N2 (5-8) 

Covariance matrix = R^ = E [ (N - p^) (N - y^) ] (5-9) 

Variance = (5-10) 

Autocorrelation function = (5-11) 



In computer simulated images, the clutter noise is 
described by either the first order or the second order 
Markov models described in Appendix A, and the PDF is 
assumed to be Gaussian. 

2 . Single Frame Processing - Nonrecursive 
Statistical Spatial Filters 

Two nonrecursive statistical linear spatial filters 
were developed to process the image signal X. 

The filter is described by its filter coefficient 

vector : 
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(5-12) 



W = [ w^ w 2 • . • Wj ] 

The coefficients can be designed by two criteria. A 
minimization of mean square error led to the spatial MMSE 
(or Wiener) filter. A maximization of output signal to noise 
ratio led to the spatial matched filter. 

The filtered output is described by: 




(5-15) 



which contains both the filtered target signal y^ and filtered 
clutter noise y^. 

The filtered statistical targets are described by the 
following statistical characteristics: 



Mean = y v 

Y S 

. 2 

Variance = a v 

Y S 



y s -W T 

w.r ss -w t 



(5-14) 

(3-15) 



The filtered deterministic targets are simply described 
by their intensities. 

y s = S-W T (5-16) 

The filtered clutter noise is described by the follow- 
ing statistical characteristics: 

1 • T 

Mean = p Y - Iw(i)u M (i) = y v *W (5-17) 

Y N 1 iN N 

Variance = o?, (5-18) 

y N 

= E[(y N -PY N ) T ( Y -VY N )] = E[W(N-y N ) T (N-y N )W T ] 
= W E[(N - y N ) T (N - u N ) ]W T 
= w R NN w t 
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3 . Multiple Frame Processing - Nonrecursive 
Statistical Spatial- Temporal Filters 

In multiple frame processing, a one dimensional 

temporal filter is added to the two dimensional spatial 

filters W. In this thesis, spatial filtering is performed 



first on several frames of images 
temporal filtering of the spatial 
nals from the same pixel location 
frames. This set of time-varying 
the vector 

Y K = ^ y l, y 2, ' * ,y K^ 
consists of target signals - 

S K = t s Kl, s K2, ' ’ ’ S KK^ 
and correlated temporal noise - 

N K = ^ n Kl , n K2, • * ’ n KK^ 



It is then followed by 
filtered images using sig- 
from several consecutive 
signals is described by 

(5-19) 

(5-20) 

(5-21) 



In the computer simulations reported later, only the first order 
Markov model is used for correlated temporal noise. The 
temporal filter is described by: 

W K t w Kl, W K2, ‘ ‘ ' W KK^ (S' 22 ) 



The output spatial - temporal filtered image is described 
by: 

^ = y k ‘ w k T C 5 - 23) 

Again, z also has its target components and clutter noise 
components . 
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The spatial- temporal filtered target signals are 
described by: 

T 

The signal intensity = Zg = ^K*^K (5*24) 

for deterministic targets. For statistical targets, they 
are described by: 

Mean = yz<, = (5-25) 

2 T 

Variance = 0Zg = ^K*^SSK*^K (5*26) 

C. TARGET DETECTION IN FOCAL PLANE PROCESSING 

Two methods of calculating threshold level are investi- 
gated . 

1) CFAR (constant false alarm rate) 

2) Bayes criteria. 




Figure 5-1 

1. CFAR 

In the CFAR approach, the threshold level = is 

selected such that Pp^ is constant. It is important in many 
practical military applications because a system will operate 
for a long time without any target being present, facing 
mostly background clutter noise. It is important to know 
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the expected false alarm rate and keep it constant. Two 
methods of setting will be presented, 

a. Gaussian PDF Case 

When the PDF of the filtered noise (denoted by 
P(u)) is Gaussian, as shown in figure 5-1, the 
threshold level can be calculated analytically as follows 



P 



FA 



QO 



/ P(u)du 
^TH 



CO 

/ P (u) du 

- CO 




P (u) du 



(5-27) 



as shown in Fig. 5-1 by the shaa.ec area. 

u denotes noise variable 
u u - E[uJ 

°n m E [fu-w u ) 2 ]. 

Since P(u) is a PDF, the first term is 1. The second term 
is the normal cumulative distribution function F(u) which is 
tabulated or can be computed using standard computer routines 



FA 



= 1 - F (!™_Ji ) 



u 



This threshold equation is: 



V TH = P FA> 



u 



(5-28) 



a u can be found from a u which is the variance of the 
filtered noise given by (5-18) for the spatial filter case 
and by (5-26) for the spatial- temporal filter case. y u is 
the mean of the filtered noise given by (5-17) for the 
spatial filter case and by (5-25) for the spatial- temporal 
filter case. 
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In many practical cases, the mean of the filtered 
noise is approximately zero. Therefore the threshold level, 



V TH - K- °u C 5 - 29 > 

is proportional to the standard deviation of the filtered 
noise for a specified PpA* For CFAR requirements, the setting 
of is equivalent to finding the constant K. 

b. Non-Gaussian PDF Case 

In real world situations, it was found that the 

PDF of clutter intensity is usually non-Gaussian. A typical 

example is shown in Fig. 5-2 which contains normalized 

histograms of real world infrared scene. 

Fig. 5-2(a) - "ApDF" of unprocessed image ("APDF” is approximated PDF) 

Fig. 5-2(b) - "APDF" of spatial - filtered image 

Fig. 5- 2(c) -"APDF" of spatial- temporal filtered image. 

The "APDF" is achieved by evaluation of each pixel intensity 
and counting all pixels with approximately the same intensity 
so as to achieve the histogram. The histogram is then normal- 
ized so that the total area under the curve becomes 1 so as 
to yield the approximated PDF. 

It can be seen that the mean of the filtered scene 
can be approximated by zero, and also that none of these "APDFs" 
are Gaussian. 

For these cases, a numerical algorithm is devel- 
oped to set the threshold for a specified Pp^, or CFAR. The 
basic assumptions are the following: 
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figure 5-2 

Typical Histograms of Clutter Intensity 
(Hand fitted curve is added to the histogram to indicate au approximate PDF) 



a) The "APDF" of filtered image is available but cannot be 
expressed in terms of an explicit analytical function. 

b) The mean of the filtered image is zero. 

c) Similar to the Gaussian PDF case, the threshold level 
is given by: 

V T „ = K*o (a is standard deviation 
u u of filtered image) 



K is to be determined by a computational procedure whose 
flow chart is shown in Fig. 5-3(a)(b). The false alarm rate 
can be shown to be indeed constant for a given K considering 
the following: 



P FA “ ( P{U)dU 
V TH 



(5-50) 



where 
to be 



p(u) is any given PDF. The threshold level is selected 



V TH = K-a 



(5-31) 



where 



7 



= / u“p(u)*du = Cp (assuming u has zero mean) 



(5-52) 



If we select C 0 so that a = /TT^ = C-, and substitute eq. 
(5-23) into eq. (5-22), the Pp^ for a given K is: 



FA 



+ oo 

f p Cu) du 
K-Co 



= constant = C. 



(5-33) 



K=constant 



Eq. (5 -25) shows that for a given K, 
the CFAR algorithm. In general, Pp^ 
can be expressed as Pp^(K) . 



FA 
is a 



is constant, 
function of 



yielding 
K and 
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Figure 5-3[a) Procedure for Setting Threshold in Cases 

of Non-Gaussian Probability Density Function 



104 



PFA(k) 
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Flow chart of the 
generation of vector 
PFA(K) 



(let next 
pixel (m,n) 
= x ;X=x/o 




vector; 
PFA(K) + 1 ; 
K=K+A 
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Figure 5- 3(b) Flow Chart for Generating PFA versus K Curve 
Used in Setting Threshold 
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The only nontrivial step in this flow chart is 
step 3 where we generate a vector (or table look up) which 
relates K to P pA [vector Pp A (K)] for a given a u and unknown 
PDF. Figure 5-4 is a typical plot of such a table look up 
for Gaussian PDF. In Fig. 5-5, the vector P p ^(K) will con- 
tain data similar to Fig. 5-4 but for the unknown PDF. To 
perform step 4, a P p ^ has to be specified and then the vector 
P p ^(K) will be searched (left to right) to find the first K 
to yield P p ^ <_ P p ^ (specified). 

2 . Bayes Criterion 

Bayes criterion is another method of threshold setting 
based on a priori knowledge of the target as well as the noise. 
The basic disadvantage of this method, compared with the CFAR 
approach, is the requirement to specify the target, background 
clutter and their statistical properties. However, if this 
information is already available, the Bayes criterion can be 
useful. Otherwise, the CFAR approach is more useful. The 
Bayes criterion is described in detail in Appendix B. 

3 . Probability of Detection 
a. Introduction 

If the probability density function (PDF) of the 
target and also the noise are known, the probability of detec- 
tion (P^) can be calculated in a straightforward manner. It 
should be emphasized that the clutter noise considered in the 
calculation of P^ is, in most cases, not the unfiltered 
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*PDF(S) of Signal Alone 



Figure 



Figure 





5-9 Probability Density Function of Signal + Noise 

(a) Target Signals 

(b) Noise 

(c) Targets and Noise 



Fig . 



5-10 




Calculation of Probability of 
Detection 
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clutter, but the filtered clutter noise. Because in the 
unfiltered condition intensities of targets are usually only 
a small fraction of the clutter intensity, the threshold 
should not be applied at this stage. Only in limited cases 
is target intensity sufficiently higher than the clutter so 
that a threshold can be applied without any predetection 
processing . 

Two methods of calculating the are implemented 
in this thesis. The first method is analytical and only 
applies to the special cases when the PDF is either Gaussian 
or can be explicitly described by mathematical expressions. 
The second method is numerical and is applicable to any case, 
regardless of whether or not the PDF can be described by well 
known functions. This method is implemented in the computer 
simulation by counting the detected targets and comparing 
their number with the number of inserted targets. 

b. Gaussian PDF Case 

If the PDF of the clutter noise is Gaussian, the 
PDF of the processed clutter by any linear filter is also 
Gaussian, as shown in Fig. 5-9(b). Target signals considered 
in this thesis can be either deterministic or statistical, as 
shown in Fig. 5-9(a). Assuming the signal and noise are 
uncorrelated, then 

PDF(s + mj = PDFf g-j *PDF is also Gaussian C5-34) 

where * denotes convolution. A typical result is shown in 
Fig. 5-9(c). The resultant statistical properties are: 
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y S + N 

„ 2 
a S+N 



y + u M = mean of signal + noise 

(5-35) 

2 ? 

Og + = variance of signal + noise 



The mean and the variance have been calculated 
in Chapter V,. section B. 

For a selected threshold , the probability 
of detection is given by the shaded area in Fig. 5-10 and 
can be determined by the following integration: 



P d = / 



V TH a S + N 



exp [ - T ( 



1 / X ' y S+N 2 



S+N 



) J dx 






/ N ^S+N, a S+N> d C S+N > ' ; ^S+N^S+N^ 5 ^ 



P 



d 



1 - 



F ( 




~ y S + N ^ 
S+N 



(5-36 ) 



x 

where F(X) = / N(0,l)dn 

— oo 



= normal cumulative distribution function. 

F CX) can be easily evaluated by using the error function 
routine available in most computer libraries because 

erf (X) = 2F(/2x)-l (5-37) 

c. Non-Gaussian PDF Case, P^ 

If the PDF of the unprocessed clutter noise is 
other than Gaussian, the PDF of the processed clutter is 
generally not Gaussian. The processed clutters PDF can be 
obtained wither analytically or by a best fit algorithm, and 
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the PDF of the processed signal plus clutter can be obtained 
as well. Denoting the PDF of clutter as P (N) and the PDF of 
signal plus clutter as P(S+N), the resulting probability of 

CO 

detection is given by P, = / P(S+N)d(S N). 

a V 
V TH 

D. SUMMARY 

Three major steps in target detection are discussed in 
this chapter. First, the probability density function (PDF) 
of both filtered image and target signals are determined. 

Next, based on the PDF, a threshold can be calculated for 
any specified probability of false alarm. Third, from the 
threshold, the probability of detection can also be calculated 
using the known PDF. 

If the PDF of filtered clutter noise is Gaussian, both 
the threshold and P^ can be theoretically determined. Other- 
wise, they must be numerically calculated using the known 
PDF. Typical results based on computer simulations will be 
presented in Chapter VI where predetection processing using 
statistical spatial filters and temporal- spatial filters for 
clutter suppression will be combined with thresholding for 
target detection. Both computer simulated images containing 
Markov model clutter noise and a real world infrared back- 
ground image will be used to test the results of focal plane 
processing . 
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VI . PROPERTIES AND PERFORMANCE OF NON ADAPTIVE 
STATISTICAL FOCAL PLANE PROCESSING - 
NONRECURSIVE FILTERS FOR CLUTTER SUPPRESSION 
FOLLOWED BY THRESHOLDING FOR TARGET DETECTION 

A. INTRODUCTION 

Focal plane processing consists of many steps. It is 
started at the optical detectors. After detection, multi- 
plexing and signal conditioning, the first two steps are 
clutter suppression and thresholding. Following these, image 
enhancement techniques may be used to improve the shapes and 
features of targets, linking procedures may be used to con- 
nect moving point targets from many frames to form tracks in 
what is sometimes known as the "track association" or "track 
assembly" processing, trajectory estimation techniques may 
be used to determine the dynamics such as velocity and 
acceleration of the targets, pattern recognition and discrim- 
ination techniques may be used to obtain more information on 
the targets. All of these additional processing techniques 
can be used to distinguish true targets from false alarms. 
However, this thesis is concerned only with the first two 
steps - background clutter suppression and thresholding. 

In previous chapters, the principles and design proced- 
ures for four focal plane processing procedures have been 
developed. 
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For background clutter suppression: 



Chapter III - Single frame focal plane processing: 
MMSE filters 
Matched filters 

Chapter IV - Multiple frame focal plane processing: 
MMSE filters 
Matched filters 

0 For initiation of target detection: 

Chapter V - Thresholding 

In this chapter, quantitative results of these focal 
plane processing procedures will be investigated and analyzed 
in more detail. The organization of results of computer 
simulated analysis and tests are described in the following 
section . 

B. ORGANIZATION OF COMPUTER SIMULATED STUDIES 

1 . Types of Focal Plane Processing Procedures 

Four focal plane processing procedures have been 
developed, as listed above. Two are spatial filters which 
use one frame of an image. Two are spatial- temporal filters 
which use several frames of images. Each of these filters 
is followed by a thresholding step. 

The important features of these filters are: 
a. For spatial filters; 

0 3x5 search box 
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Estimation point is in the center of the box 
in most cases. It is moved to the edge or 
corner of the search box when the estimation 
point is along the edges or corners of an image. 
0 Point targets are considered in almost all cases 
b. For temporal filters: 

0 One pixel from three or five successive frames. 

0 Estimation point is in the middle of successive 
frames. Example: second frame for a 5 frame 

temporal filter, third frame for a 5 frame 
temporal filter. 

0 Point targets are considered which move at 
least one pixel from one frame to the next. 

2 . Types of Images Used in Computer Simulation 
Two classes of images are used, 
a. Computer Generated Images 

Unprocessed images are generated by computer 
simulation and contain three components: 

(1) Targets - Mostly point targets, described 
by their locations and intensity. 

(2) Uncorrelated random noise is described by 

2 

its mean, which is usually zero, and variance . 

(3) Correlated random noise - Two types of Markov 
processes are used to describe correlated random noise. The 
detailed theoretical model analysis is presented in Appendix A 
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First order: Q(m,n) = a?, N P^Py^ 

Second order: Q(m,n) = p^ m ^ p| n ^ cos C8 H *m) cos (By • n) 

It can be seen that three parameters are needed to 

describe the first order Markov model. 

2 

a CN = var i- ance °f correlated noise 
Ppj = horizontal correlation factor 
Py = vertical correlation factor. 

For the second order Markov model, two additional parameters 
are needed. 



Bpj = horizontal spatial frequency 
By = vertical spatial frequency 



b. Real World Infrared Image 

Twenty successive frames of a real world infrared 
image in the medium wavelength infrared spectral band are 
used in the test of sequential spatial- temporal filters. 

The target intensity and the unprocessed 
image are combined together and described by another very 
important parameter. 



Input signal to noise ratio = SNRI 



3 . Properties and Performances Evaluated 

The following properties and performance figures of 
merit are evaluated. 
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a. For Clutter Suppression: 

(1) Filter coefficients as functions of input 
signal to noise ratio and correlation factors. 

(2) Processing gain as functions of input signal 
to noise ratio and correlation factors. 

(3) Spatial side-lobes after the application of 
spatial filters as functions of input signal to noise ratio. 

b. For Threshold Detection: 

(1) Threshold levels as functions of input signal 
to noise ratio, using the probability of false alarm (Pp^) 

as a parameter. 

(2) Probability of detection as functions of 
input signal to noise ratio and probability of false alarm. 

4 . Sensitivity Analysis 

All nonrecursive filters developed in this thesis 
use a common set of mathematical operations which is the 
summation of weighted and delayed signals. The differences 
among four nonadaptive filters developed in Chapters III and 
IV and adaptive filters to be developed in Chapters VII and 
VIII are in the criteria and design procedures of their 
filter coefficients. Therefore, there exists a real interest 
in using new sample analog tapped delay line devices which 
can efficiently carry out the nonrecursive filter operations 
at high speed because it is a pipelined and parallel processor. 
However, the disadvantage of sampled analog tapped delay line 
implementation is its lower accuracy. Therefore, a special 
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investigation is carried out by truncating and rounding-off 
the filter coefficients to simulate the lower accuracy of 
sampled analog implementation and investigating their 
effects on filter performance. 

5 . Suboptimal Filter Investigation 

In practical applications the noise parameters that 
were used to design the filter coefficients might be in- 
accurate, resulting in a "suboptimal" filter. The suboptimal 
filter performance is investigated in this thesis by monitor- 
ing the PG of a given "suboptimal" filter while the background 
noise correlation factor (p^ = Py) is changed from 0 to 1. 

C. SINGLE FRAME PROCESSING I - NONRECURSIVE SPATIAL 

MMSE FILTER FOLLOWED BY THRESHOLDING 

1 . Introduction 

An isotropic second order Markov process model is 
used to simulate the correlated spatial noise in this section. 
The spatial frequencies of the correlated noise are 
$h = By = ^ = 0.267. The horizontal and vertical correlation 
factors are varied. The variance of the correlated noise is 
1.0. The variance of the uncorrelated noise is 0.01. The 
target intensity is also varied and is described, together 
with the spatial noise, by the parameter "input signal to 
noise ratio" (SNRI). 

2 . Filter Coefficients 

Since the unprocessed images contain only isotropic 
random noises and point targets, the filter coefficients are 
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centrally symmetrical. Of the nine filter coefficients, 
there are only three distinct values: center coefficient W^, 

corner coefficient = W- = = Wg , and edge coefficient 

^2 = ^4 = ^6 = ^8* They are presented in Fig. 6-l(a) and (b) 
as functions of input signal to noise ratio and for five 
different correlation factors, = Py = p = 0.999, 0.99, 

0.9, 0.8 and 0.4. It should be pointed out that p = 0 cor- 
responds to the uncorrelated random noise case and p = 1 
corresponds to the deterministic case. 

The following observations of their behaviors are 

made . 



a. Symmetry 

Filter coefficients have a central symmetry as 
described above. 

b. Coefficient Values 

All coefficient values are less than unity. They 
are easy to implement in hardware if sampled analog image 
processing is used by simply using a resistor potentiometer. 

The sum of nine coefficients is very close to 
zero. This implies that this type of spatial filter will 
usually suppress the spatial mean of background noise. It 
can be considered as a "whitening" spatial filter. 

c. Dependence on Correlation 

As the correlation factor approaches unity, the 
filter coefficients change less and less as p varies and are 
also independent of input signal to noise ratio. In the limit 
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Miter coef f i c ient s conf igurat i on . 



of p = 1 , the filter coefficients approach the values of a 
spatial Laplacian mask in those cases when 3^ and By are 
small C< 0.26). 

The Laplacian 0.25 -0.5 0.25 

mask for point 

target : -0.5 1 -0.5 

0.25 -0.5 0.25 



For high values of spatial frequencies, 3^ and By (>0.26), 
the filter coefficients are different than the Laplacian mask. 

In the other direction, when p approaches zero 
all eight filter coefficients other than the middle coeffi- 
cient are decreased. In the limit of p = 0, only the middle 
coefficient is non zero. This indicates that in the extreme 
case of uncorrelated spatial noise, there is no help provided 
by using neighboring signals to estimate a point target and 
to suppress the background noise because the noise is not 
correlated . 



d. Dependence on Input Signal to Noise Ratio 

Filter coefficient values become smaller as the 
input SNR is decreased. It can be understood because the MMSE 
filter is designed to give the optimal estimate of the target 
signal. When the noise is high, the filter tends to attenuate 
the noise by smaller coefficient values and to yield the 
estimation of the small target value. 
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Processing Gain 



3 . 

Processing gain of the spatial MMSE filter is pre- 
sented in Fig. 6-2(a) as functions of input signal to noise 
ratio for different correlation factors and in Fig. 6-2 (b) 
as functions of correlation factors. The following obser- 
vations are made. 

a. Dependence on Input Signal to Noise Ratio 

The processing gain of spatial MMSE filter for 
point target estimation is independent of input SNR. 

b. Dependence on Correlation 

The processing gain is strongly dependent on the 
correlation factor. The higher the correlation, the larger 
processing gain can be achieved, but this will be affected 
by the intensity of uncorrelated noise components also in the 
image. For example, in this particular case, the PG of approx- 
imately 45 db is achieved for p = 0.999, as shown in Fig. 6-2(a). ■ 
If the uncorrelated spatial noise is excluded, which is the 
case in Fig. 6 - 2 ( b ) when p = 1, the PG can be as high as 60 db . 
Unfortunately, uncorrelated noise always exists in real world 
images due to background shot noise, electronic noise, etc. 

As approaches zero, the PG decreases. In the 
limit p = 0, PG = 1 which confirms the discussion in the last 
section on filter coefficients. 

4 . Threshold Level and Spatial Side Lobe Level 

Threshold level and spatial side lobe level are 

related because if the threshold is set below the side lobe 
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Figure 6-2. Processing Gain of Nonrecursive Spatial Filter as 
a Function of X a ) Input Signal to Noise Ratio, 
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level, a false alarm will result. The spatial side lobe 
level is determined by multiplying the highest positive 
filter coefficient with, the target level. 

In Fig. 6-30), (b), threshold levels are plotted as 
functions of input SNR for several values of probability of 
false alarms from 10 ^ to 10 ^ and a given correlation factor. 
Curves are presented for p = 0.99, p = 0.9, p = 0.8, and 
P = 0.4. Two observations are made: 
a. Dependence on Correlation 

As P is increased, the threshold value can be 
lowered to achieve the same Pp^. This is because the pro- 
cessing gain of the spatial MMSE filter is better for higher 



P. 



b. Dependence on Pp ^ 

As Pp^ is reduced, threshold value must be in- 
creased . 



The spatial side lobe levels are presented in 

Fig. 6-5(a),(b) as well as functions of input signal to noise 

ratio. By comparing the threshold level plots and side lobe 

level plots corresponding to the case of p = 0.9 in Fig. 

6-3(a), it can be seen that false alarms due to side lobes 

will occur for input SNR approximately above -4 db for 

- 2 - 1 

threshold levels in the range 10 “ < Pp^ < 10 
5 . Probability of Detection 

A set of figures (Fig. 6-4) is used to present the 
probability of detection, P^, as functions of input SNR for 
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Figure 6-3(b) 
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several values of Pp^ from 10 ^ to 10 Each figure is for 

a selected correlation factor, = Py. 

It should be emphasized that the two figures of merit 
of PG and P^ must be used together in the evaluation of 
effectiveness of focal plane processing. Otherwise, mislead- 
ing evaluation could result as demonstrated by the following 
example . 

Consider a case of p = 0.9; the processing gain of 
a 5x3 spatial MMSE filter for point target estimation is 20 
db. It could be interpreted that any target of input SNR 
better than -20 db, say -12 db , can be detected. However, 
the P^ plot in Fig. 6-4 reveals that for Pp^ = 10~^, the 
probability of detection of a weak target of -12 db is very 
low. P^ is less than 0.9 even for a stronger target of -7 db . 
Therefore, the capability of a focal plane processing process 
consisting of a spatial filter for clutter suppression (eval- 
uated by PG) and thresholding (evaluated by Pj) must be cor- 
rectly determined by using both PG and P . 

6 . Summary 

Some highlights of the first type of focal plane pro- 
cessing which consists of spatial MMSE filtering and thresh- 
olding are; 

a. The processing gain is strongly affected by 
correlation factor, relative intensity of white noise, types 
of targets, search box size and location of estimation pixel 
[Ref. 14J . 
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b. Input SNR does not affect the processing gain 
but will influence the probability of detection. 

D. SINGLE FRAME PROCESSING II - NONRECURSIVE SPATIAL 

MATCHED FILTER FOLLOWED BY THRESHOLDING 

1 . Introduction 

The same unprocessed image used in the study of 
nonrecursive spatial MMSE filters is used in this section. 

There will be a great deal of similarities in the organiza- 
tion of presenting the results also. 

2 . Filter Coefficients 

Figure 6- 5 (a) , (b) , (c) presents the filter coefficients 
as a function of input SNR for several isotropic correlation 
factors. Since point targets are being considered and only 
isotropic noises exist in the unprocessed image, the filter 
coefficients have again a central symmetry. Only three of 
the nine filter coefficients need to be presented - W(l), 

W(2) and W(5) . 

It can be seen that their values are not only very 
large but also have strong dependence on the input SNR. Both 
facts are inconvenient in hardware implementation. (It should 
be remembered that in spatial MMSE filters, the filter coeffi- 
cients are all less than one and are independent of input SNR.) 
These difficulties with the matched filter can be removed by 
a proper choice of the constant C in the solution of the 
coefficient equation C^-4.28) to normalize the filter coeffi- 
cient for the estimation pixel to unity. Figure 6-6(a),(b) 
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Fig. 6-6(a) Matched Filter Coefficients as a Function of Input SNR; P,. = P v = Parameter 
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filter Coefficients Configurations 



shows the dependence of filter coefficients after the normal- 
ization. 

Other observations made for the MMSE filters also 
apply for the matched filters. 

3 . Processing Gain 

The processing gain of the spatial matched filter for 
point target estimation is presented in Fig. 6-7 as a function 
of input signal to noise ratio for different isotropic corre- 
lation factors, pj_j = py 

Again, the general observations obtained for the MMSE 
filters are also valid in the matched filter cases. In fact, 
by comparing Fig. 6-7 with Fig. 6-2(a), it is discovered that 
although the filter coefficients for these two filters are 
quite different, their processing gains are the same. This 
result can be justified by using the following reasoning. The 
MMSE filter is designed to suppress the noise and to yield the 
best estimate of the signal. The matched filter is designed 
not only to suppress the noise but also to enhance the signal. 
This second feature of the matched filter results usually in 
signal distortion in contrast to the MMSE filter which is 
designed to yield the best estimate of the signal. When the 
concerned signal is a point target, as in the above case, the 
normalized matched filter cannot affect the signal and there- 
fore both matched filter and MMSE filter will only suppress 
the noise to the same extent, resulting in the same PG . 
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4 . Threshold Level and Side-Lobe Level 

Threshold levels are plotted in Fig. 6-7(a),(b) as 
functions of input SNR for several values of Pp^ from 10’ 1 
to 10 and a selected correlation factor. Several figures 
are presented for different p's. 

Side lobe levels are also shown in Fig. 6-7(a),(b) 
as functions of input SNR. 

Again, there are many similarities between the spatial 
matched filters and the spatial MMSE filters presented in the 
last section. It should be noted that the threshold level of 
the normalized matched filter does not vary with SNRI as it 
does for the MMSE filter. 

5 . Probability of Detection 

The probability of detection is presented in Fig. 6-3 
as functions of input SNR for several values of Pp^ from 10 ^ 
to 10 ■*. Each figure is for a selected correlation factor. 

6 . Summary and Comparison with MMSE Filters 

Comparing spatial matched filter results of this sec- 
tion and spatial MMSE filter results in the previous section, 
many similarities are observed. The only difference is in 
the filter coefficients. It must be pointed out that these 
observations are true only for point targets. For other 
targets, MMSE filters and matched filters will not be as 
similar as in the case of point targets. It is believed that 
for extended targets, the matched filters will offer more 
superior processing gain and probability of detection because 
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Fig. 6-7(l>) Threshold Level (V^) and Spatial Side Lobe Level (SLL) as a Function of 
Input Signal to Noise Ratio in db (SNRIfdb]); PFA = Parameter 
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the matched filter is designed to enhance the target as well 
as suppressing the noise, as mentioned in VI. D. 3. 

E. SUBOPTIMAL SPATIAL FILTERS 

1 . Filter Coefficients Accuracy 

When filter coefficients are inaccurate due to round- 
off error or other implementation constraints, the processing 
gain will deteriorate accordingly. Fig. 6-9 presents the 
deterioration in PG as a result of artificial coefficients 
round-off at a certain location after the decimal point. It 
is seen that round-off at the fourth place after the decimal 
point hardly affects the PG . 

2 . Suboptimal Filter Performances 

When the correlation factor of the noise is expected 
to vary within the image, it is necessary to determine which 
is the best "suboptimal" filter. Figure 6-10 shows the 
resulting PG when a filter designed for a correlation factor 
of p = 0.999 is facing a correlation factor that varies from 
0 to 1. Comparison of Fig. 6-10(a.l) with Fig. 6-9(a.l) 
reveals that the resulting deterioration in PG in various p's 
is from 0 to about 5 db . Figures 6- 10 (a . 2) , (a . 5) , (a . 4) and 
(a. 5) are similar to Fig. 6-10(a.l) with filters designed for 
p=0.99, 0.9, 0.8 and 0.4 respectively. 

The strongest deterioration in PG is noticed in the 
suboptimal filter related to p=0.4 (Fig. 6-10(a.5)) and the 
least deterioration is noticed in the filter related to 
P = 0.999 (Fig. 6- 10 (a. 1) ) . 
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big. 6- 10 (a) Processing Gain [db] as a function of p - of suboptimal filters. 




;6 



-n 

r- 



l n 



uo 

cn 



o 


rr 


CM 






LO 


r- 




r^ 


r^ 


O 


o 


LO 




to 


CM 








141 



Fig. 6-10(b) Processing gam |<Jb] ns a function of p suboptimal filter. 



Conclusion: The best suboptimal filter in this case 

was found to be the one designed for p= 0.999. Similar sub- 
optimal studies should be performed to select the best "sub- 
optimal" filter whenever noise statistics are expected to 
vary and when filter coefficients cannot be adapted on line. 

F. MULTIPLE FRAME PROCESSING - NONRECURSIVE SEQUENTIAL 

SPATIAL- TEMPORAL FILTER FOLLOWED BY THRESHOLDING 

1 . Introduction 

To evaluate three dimensional spatial- temporal focal 
plane processing, successive frames of images are needed. 

In Chapter IV, theoretical models based on Markov processes 
were used for correlated noises. In the spatial domain, 
both first order and second order Markov model correlated 
spatial noises are considered. Only first order Markov model 
was considered for the temporal correlated noise. 

In this section, instead of using those computer 
generated images containing Markov model correlated noise, 
one set of real world infrared images of 32 x 32 pixels with 
constant drift of 10% of a pixel between frames is used to 
evaluate the performance of the sequential spatial- temporal 
filter followed by a thresholding. Except that it was orig- 
inally taken by an ERIM 0 background measurement and subse- 
quently modified by another group to include drift and several 
other modifications, no detailed information was given because 

°ERIM is the Environmental Research Institute of Michigan 
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of the classified nature of the program which developed this 
scene as a standard test for focal plane processing algorithms. 

In this section, statistical spatial filters are 
first designed for each frame to suppress clutter. The sta- 
tistical spatial correlation properties of every frame of 
images are calculated separately and used in the filter 
design. Then, several frames of these spatially-filtered 
images are used together for the calculation of their temporal 
correlation properties which are used to design a temporal 
filter which will be applied to several frames of spatial- 
filtered images to give one frame of spatial- temporal filtered 
image. A threshold level is then determined from the PDF of 
the intensity of the filtered image and applied to initiate 
the target detection. These procedures were developed in 
another thesis by Captain D. Hilmers [30]. 

2 . Spatial Filters and Their Performance 

Figure 6-11 is a nine grav-level* computer overlay 
print of one frame of the real world infrared image. Figure 
6-12 is a plot of its probability density function. The mean 
is found to be 4.03 and the variance is .45. Six moving 
targets of intensity 0.001, shown in Fig. 6-13, are added to 
this frame of image. This corresponds to an input signal to 
noise ratio of -72.3 db . The targets move down 1 pixel for 
each frame. 

*Nine gray level: "blank" -,+, V, M, 0, 9, 9, E, M 
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Figure 6-11 9 Gray-level Computer Print of One 

Frame of 32 by 52 Pixels "Mill Creek" 
Infrared Image. 
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1 lie Histogram ol the Intensity of the Original Clutter linage Normalized so that 
the Area under the Curve is 1.0. [A hand fitted curve is added to indicate the 
shape of the PDF] 
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2.20 3.27, • , , ^-35 , m 5.42 . (v. 50 

Relative Background Clutter Intensity 

Figure 0-12 Histogram of Background Clutter intensity of "Mill Creek" Image. 
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Figure 6-13 Six Moving Point Targets Added 

to the First Frame of "Mill Creek" 
Background Image. 
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At this low input SNR and for point targets there 
is no difference between the MMSE filter and the matched 
filter. A 3x3 spatial filter is designed and its filter 
coefficients are shown in Fig. 6-14. This spatial filter 
is then applied throughout this 32 x 32 image. The spatial- 
filtered image is shown in Fig. 6-14 with its new PDF shown 
in Fig. 6-15. It can be seen that the new PDF is somewhat 
symmetrical around zero intensity, indicating a zero mean. 

3 . Temporal Filter and Its Performance 

Five consecutive spatial filter frames of images are 
filtered by a temporal matched filter whose coefficients are 
shown in Fig. 6-16. The resulting temporal filtered image 
is shown in Fig. 6-16 with its PDF shown in Fig. 6-17. The 
processing gains of both the spatial and temporal filters are 
also listed in Figs. 6-14 and 6-16 respectively. A threshold 
is then set to yield a Pp^ of 10 0 and the resulting image is 
shown in Fig. 6-13. It should be recognized that the third 
frame is the "estimation frame" of this five frame temporal 
filter. Consequently, the targets in Fig. 6-18 appear 3 pixels 
down relative to their location in Fig. 6-13. 

4 . Conclus ions 

The spatial- temporal filtering method is a very useful 
technique for suppression of correlated clutter and detection 
of low signal level target. In the general case when K two 
dimensional frames of images are given, an optimal spatial- 
temporal filter can be designed for a given shaped target to 
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Spatial Filter Coefficients with 
Asymmetry Due to Anisotropy of 
Statistical Properties of the 
Real IR Image. 
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Figure 6-14 Nine Gray-level Computer Print of 
Spatially Filtered First Frame of 
"Mill Creek" Image 
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Jhe Not inai i zed Histogram of Spatially Filtered Image [A hand 
fitted curve is added to indicate the shape of the PDFJ 
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The Spatially Temporally Filtered Image 
and the Temporal Filter Coefficients 

Temporal Filter Coefficients 
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Figure 6-16 Nine Gray-level Computer Print of 
Spatially and Temporally Filtered 
Third Frame "Mill Creek" Image Including 
Background Clutter and Point Targets. 
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The Normalized Histogram of Spatially Temporally Filtered Imag 
1A hand fitted curve is added to indicate the shape of the PDF 
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Figure 6-17 Histogram of Spatially and Temporally Filtered Third Fr 
of "Mill Creek" Image Including Background Clutter and 
Point Targets. 
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Figure 6-18 



Nine Gray-level Computer Print of the 
Third Frame of "Mill Creek" Image after 
Spatial -Temporal Filtering and Thresholding. 
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achieve a certain processing goal. These goals and the 
resulting filter design method have to be evaluated for the 
particular application. In this chapter, a spatial temporal 
filter for the enhancement of moving point targets has been 
designed and tested using real IR clutter images. The 
results indicated that this type of filter can yield a 
processing gain as high as 60-70 db. 
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VII . ADAPTIVE SPATIAL FILTER I - 

STATIONARY CLUTTER STATISTICS 



A. INTRODUCTION 

In previous chapters, two types of filters were discussed: 
the two dimensional and three dimensional MMSE filters and 
matched filters, along with their design and implementation 
schemes. The performance of these filters was analyzed and 
presented. They both offer significant suppression of clutter 
noise when properly designed. The first problem we face when 
trying to design the filter is the unknown statistical proper- 
ties of the noise. The second and more serious problem is 
the nonuniformities of the statistical properties of the noise. 
In other words, the noise is nonstationary. 

In this chapter, a two dimensional adaptive filtering 
technique will be developed for cases where statistical 
properties of the noise are unknown but does not vary (sta- 
tionary statistics). Extension to the nonstationary statis- 
tical properties of the image will be presented in Chapter 
VIII. An adaptive algorithm will be developed first for the 
Wiener filter, followed by a suggested modification to the 
matched filter. 

B. ADAPTIVE WIENER FILTER 

The requirement for this adaptive algorithm is to change 
the filter coefficients to minimize the mean square error as 
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fast as possible based on the clutter noise characteristic. 
The adaptation is achieved by utilizing the steepest descent 
algorithm together with a proper selection of required refer- 
ence performance. Steepest descent is a well established 
technique for one dimensional adaptive filters. Extension 
of this technique to two dimensional adaptive filters will 
be presented in the following sections. An implementation 
scheme will be suggested along with simulation results and 
performance analysis. 

It should be pointed out that there is another important 
difference from the previously developed one dimensional 
filters. Most one dimensional adaptive filters are designed 
to operate in a white noise environment with correlated 
target samples. In this study, instead of white noise, we 
are facing a highly correlated noise. On the other hand, 
targets are often uncorrelated. For example, point targets 
are uncorrelated in both x and y variables. Horizontal line 
targets are uncorrelated in the vertical direction. 

There are two approaches to suppress the correlated 
noise . 

1. The first approach is to treat the correlated noise 
as if it were a "signal" and try to predict (or estimate) its 
value in each pixel. Reference 34 demonstrates this basic 
concept utilizing an adaptive filtering technique. Refer- 
ences 13 and 35 demonstrated again this concept by utilizing 
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a two dimensional Kalman filter rather than an adaptive 
filter. 



2. The second approach is to treat the correlated noise 
as "noise," and to design an optimal noise suppression tech- 
nique based on the statistical noise and target character- 
istics . 

In this work, the second approach is taken by direct 
adaptation of the filter coefficients so as to minimize the 
mean square error between the target and the target estimate. 
This approach was selected since it seemed to be easier for 
practical implementation. 

C. ADAPTIVE ALGORITHM BASED ON STEEPEST DESCENT 

Steepest descent is an effective adaptation technique 
and was shown to be superior [21] to other adaptive methods 
when a reference input can be specified. 

The analysis of this method is similar to the analysis 
of nonrecursive spatial filters presented in Chapter III. 

The input signal vector, Xj , is defined as 



The input signals are assumed to appear simultaneously 
on all input lines at discrete spatial locations indexed by 
the subscript j ; j is the pixel number in the image as was 




(7-1) 



The filter coefficients are 



T 

W = [w^ ,w . . . w^ . . . w_ ] 



(7-2) 
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shown before in Fig. 4-2 where j = NMAX (n-1) +m; i is the 
index describing the location of the input samples in the 
search box as in Fig. 4-2 where i = MSB * ( n - 1) + m. 

The input signals are multiplied by the nonrecursive 
filter coefficients w^,w 7 ...Wj and then summed together as 
shown in Fig. 7-1. 



Input 

signals 




Filter 

coefficients 



Figure 7-1 



The general filtering configuration is given in Fig. 7-2. 

In the following, an adaptive method will be presented to 
adapt the filter coefficients W-^...W n to achieve an optimal 
estimate of a known signal. It will be shown that at steady 
state the filter coefficient will approach the Wiener filter 
coefficients. A review of the adaptive filtering formulation 
based on references 21, 22 will be presented in Appendix C. 

The algorithm derivation is based on the following model. 
At a given pixel, the output of the filter is 

y, = X. T W = W T X. 

3 3 J 



03 ) 




Input Vector 
C 

Adaptive 
Fi 1 ter 

Coefficient s 
W 



Figure 7-2 Configuration of Nonadaptive Two 
Dimensional Nonrecursive Spatial 
F ilter . 
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Let us select a desirable reference performance, dj , 
which we seek to achieve at the steady state. For example, 
in the case of point target detection ( or estimation) , we 
seek to achieve the best estimate of a point target while 
suppressing the background clutter noise as much as possible. 
Similarly, in case of line target detection (or estimation) , 
we seek to achieve the best estimate of a given line while 
suppressing the clutter noise. 

Let us define the error at pixel location j to be 



e. = d . - y . = d . - W 1 X . 
J J J J 



( 7 - 4 ) 



It is the purpose of the adaptive process to adjust 
the weighting coefficients of the linear adaptive filter 
to minimize the mean square error e . . 

In this chapter we will make the assumption that the 
clutter noise is stationary. In the following chapter, a 
method will be suggested to deal with nonstationary clutters, 
where the statistical properties can be assumed to be 
stationary only within a small area. The detailed adaptation 
method and derivation are explained in Appendix C. 
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D. TWO DIMENSIONAL ADAPTIVE FILTERING CONSIDERATIONS 



1 . Introduction 

In this section, "non-optimal" spatial filtering 
will be described to introduce the need for adaptive filter 
ing . For this purpose, a composite simulated image is used 
which consists of four subimages simulated by computer and 
having four different statistical correlated noise charac- 
teristics as shown in Fig. 7-3. 



Upper left quadrant :- 


2nd order Markov model 
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Upper right quadrant: - 
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Lower left quadrant :- 


2nd order Markov model 
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Lower right quadrant:- 


1st order Markov process 




P H = 0.99 




Py = 0.30 
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1st 




Figure 7-5 Three Dimensional Plot of Computer Generated 
Four Quadrants "Nonstationary" Clutter Image 
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Based on the statistical properties of each quadrant, optimal 
spatial MMSE filter can be designed as follows: 
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right 
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.499 1.0 .499 

.245 -.487 .243 

However, Fig. 7-4 shows the filtered image by using only one 
spatial MMSE filter optimally designed for the upper right 
quadrant but applied throughout the whole image. It should 
be pointed out that a threshold was applied after the spatial 
filtering to obtain Fig. 7-4, It can be seen that excellent 
"clutter noise" suppression was accomplished in the optimal 
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Filtered and Thresholded Image 
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*igure 7-4 Computer Print of Nonadaptively Filtered and Thresholded 
Four Quadrant "Nonstationary" Clutter Image. 
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filter quadrant. Considerable amount of clutter noise 
remained in the other three "suboptimal" filtered quadrants. 

It would be desirable if an adaptive algorithm can be devel- 
oped so that the filter coefficients will be adaptively 
changed to new optimal values when the filter is moved into 
areas having different statistics. 

An example of such an adaptive filter configuration 
is shown in Fig. 7-5. Similar to all other spatial filters 
developed in this thesis, it uses a 3x3 search box. Signals 
of all nine pixels in the box are weighted and summed to 
give an estimate of the target signal at the central pixel. 
However, there are two major differences. First, the weight- 
ing coefficients are adaptively changed by a feedback proced- 
ure as shown in Fig. 7-5. At the beginning, nonrecursive 
spatial filtering is carried out as described previously. 

The processed target signal, y^ , is now compared with the 
selected desirable reference target, dj , by subtraction and 
produces an error signal, e.. This is then multiplied by the 
adaptive loop gain constant, p , and then by the nine corres- 
ponding unprocessed signals, X, to process a correction 
weighting coefficient, AW^ , for its corresponding coefficient, 
W. . The new filter coefficients, + AW^ , are used to perform 
a new step of spatial filtering and generate an additional 
set of corrections to the filter coefficient. Then the esti- 
mation pixel is moved one pixel, a new set of unprocessed 
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Adaptive Spatial Filter 
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Box 




Unprocessed 
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Figure 7-5 Configuration of Adaptive Two Dimensional 
Nonrecursive Spatial Filter. 
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signals, X, is read out and the adaptive filter process is 
repeated again. This adaptive process is repeatedly carried 
out according to the process described in Appendix C and 
based on the LMS adaptive algorithm, equation (C-42), until 
the steady state optimal solution is reached. The second 
difference of this adaptive spatial filter from the non- 
adaptive spatial filter is in the scanning pattern of the 
estimation pixel. It is not the conventional raster that 
the estimation pixel is moved back after the completion of 
each row and is always scanned in the same direction. Instead, 
after each row, the scan direction is reversed as shown in 
Fig. 7-5. The reason for this successive reversal of scan 
direction from one row to the next is to insure that the 
change of statistical characteristics of the clutter noise 
is minimized because in this new scan motion the estimation 
pixel is moved only by one pixel. 

More details on the implementation of adaptive spatial 
filters will be developed in the next section. 

2 . Detailed Implementation of Adaptive Spatial Filter 

Several aspects of adaptive spatial filter implemen- 
tation are described in more detail in this section: 

a. Image scanning pattern 

b. Termination of filter coefficient adaptation 

c. Derivation of desired reference d. 

d. Selection of adaptive loop gain y 
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Required number of independent pixels to achieve 
adaptation 



e . 



a. The Image Scanning Pattern 

Several image scanning patterns were tested. 
Figure 7-6 presents two possible scanning alternatives. 
Comparing the scanning pattern of Fig. 7-6(a) to that of 
7-6(b), it was found that using the pattern of Fig. 7-6(b) 
to adapt the filter did not yield the optimal filter coeffi- 
cients. This scanning pattern introduced an artificial 
transition at the end of each row and therefore affected the 
adaptive process. On the other hand, the pattern of Fig. 
7-6(a) offered a smooth adaptation without these artificial 
transitions and the resulting adapted filter coefficients 
were found to be the optimal. 

b. Algorithm for Termination of Filter Adaptation 
Several algorithms were tested for termination 

of filter adaptation. The selected algorithm is based on 
the following steps and its flow chart is presented in Fig. 
7-7. 

Cl) Test the MSE at each n* adaptation step. 

(2) Keep track of the lowest MSE = MSE m ^ n (0) 

and its corresponding filter coefficients. 

(_3) Test every N* (N>>n) adaptation step whether 

the new MSE ■ = MSE . Cl) is smaller than MSE - fO) by a 

mm mm mm^- 

*n,N are constants that can be selected either by simu- 
lation or arbitrarily, as recommended later on. 
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Image Scanning Patterns 
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Figure 7-6 Two Scanning Patterns for Performing Adaptive Spatial Fi 



F = 1 + F L %] / 1 00 




Figure 7-7. Software Algorithm to Terminate Filter Adaptation 
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certain factor F (say, 51). If smaller, then continue, else 
terminate adaptation. Other algorithms based on monitoring 
the filter coefficients error or continuous tracking of MSE 
(short term observation) did not yield satisfactory results. 
The suggested algorithm was tested in various simulated images 
and it functioned well whenever the parameters F, N, n were 
selected properly. When absolutely no a priori knowledge of 
noise statistics is given, these parameters can be selected 
for the worst case: for example, n = 36; N = 360; F = 0.1%. 

When some statistical data is available, these parameters 
can be optimized to hield fastest indication of adaptation 
termination. The information about adaptation termination 
is necessary mainly for the following two purposes: (1) an 

attempt to filter the image before reaching the steady state 
solution will yield poor noise suppression; (2) in Chapter 
VIII a particular application of the adaptive filtering 
technique will be presented by dividing the image into many 
subimages and adapting a different filter to each subimage. 
Under such circumstances and when only one central processor 
is available, it is important to detect the adaptation ter- 
mination in each suhimage before moving to the next subimage, 
c. The Derivation of the Designed Response, dj 

To derive the desired response, d^ , let us return 
to the LMS adaptation algorithm presented in eq . (C-42): 
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w. + 2]i £ . • X. 
J J J 



or 



'Vi 



w j+1 - V. ■ 2pe j -X j - iW. 



(7-5) 



where AW^ is the filter coefficients correction vector at 
adaption step j. Substitution of the error equation (02) 
into ('7-5) will yield, 

AWj = 2y (dj - W T -X^)-Xj or 

AW. 

A^W = = [d. - W 1 -X j )-X j (7-6) 

Fig. 7-8 presents typical 5x5 spatial filter and the corres- 
ponding input samples for various target types. 

Referring back to eq. (7- 6), let us expand the 
matrix for the simple point target example of Fig. 7-8(a) 

A*w 1 = x 1 {[x 1 w 1 + . . (x 5 + s)w 5 + ... x 9 w g ] - s} 

A * W S = C x 5 +S ){[ x i w i + • • • (x 5 + S)w 5 + . . . x 9 w g ] - s} 

; (7 

A*w g = x 9 ([x 1 w 1 + ... (x 5 + s)w 5 + . ..x g w 9 ] - s} 



where s is used twice on this configuration - once as 
desired response dy and once as an injected signal - 
artificially to sample Xr. 

As mentioned earlier, the adaptive process 
as a low pass filter smoothing out equation (7-7 ) to 
its expected value: 



the 

injected 

acts 

yield 
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POINT TARGET AT CENTER PIXEL 
OF SEARCH BOX " 7 -8(a) 



POINT TARGET NOT IN 
SEARCH BOX CENTER - 7-8 (b) 






AW Aw^ Aw 9 Aw 3 Aw 4 ,A^ A^ Aw., Av^ Aw g t 



Figure 7-8 
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E [ ] E [^2 X ]_] W l + ' • ^[X^Xg]Wg + • • E [x 1 X g] Wg 

E [ A’Vj.] - E [x^x^] + E[x^x^]Wg + E[x^Xg]Wg 

+ E[s-s]w 5 - E [ (x 5 + s) -s] 

* ( 7 ' 8 ) 

E [ A W g ] = E[XgX 1 ]W 1 + E[XgX 5 ]W 5 + ..E[XgXg]Wg 

Assuming the signal s and the noise samples X are uncorre- 
lated. At steady state E [Aw* ^ ] i=l , 2 , . . . 9 will be zero and 
eq. (7- 8) can be rewritten in the vector form, 

f R NN + R SS^ w = p (7-9) 

Assuming the input samples X consist only of noise and the 
target is summed artificially to X. 

R NN = C0V [ noise ] 

Rg 2 = cov [target] 

P = cross correlation vector given by eq. (C-5). 

Eq. (7- 9) is indeed the Wiener filter equation for a point 
target. We can proceed in a similar way to construct d^ for 
other simple type targets. 

Let us summarize the steps of deriving the d^ in 
the spatial adaptive filter configuration: 

1. The target shape should be drawn with respect 
to the search box dimension. 

2. All the X vector locations through which the 
target will pass will be "marked" to establish a reference 
vector, D = [d^d ? ...dj]. All the "marked" locations, d^, 
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will be set to s = signal desired level, and the remaining 
location will be set to zero according to the shape of the 
target of interest. 

3. While adapting the filter coefficients, it 
is assumed that relatively few pixels contain the target. 

An even more reasonable assumption can be made for practical 
imaging systems, considering the total absence of the external 
target during the initial adaptation phase. Making this 
assumption will require a preliminary "learning" period prior 
to utilizing the filter. Under this assumption the referenced 
vector D may be summed together with the input sample vector, 
X. The resulting outcome will be applied to the filter to 
yield its output, y^ . 

4. The output, y. , will be subtracted from the 
referenced signal d^ = s = signal level to yield the error 
signal . 



5. The error signal multiplied by the loop gain, 
2y , will be thus multiplied by each of the input vector (X + D) 
components to yield the correction vector, Aw. 

6. At steady state after many adaptation steps, 
E[aW] will approach zero and the Wiener filter solution will 
be reached. 



d. The Corresponding "Adaptive" Matched Filter 

It was mentioned earlier in this work that the 
selected major goal is to detect point targets. It is inter- 
esting to note the correspondence between the adaptive Wiener 
filter for point target enhancement and its counterpart - the 
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matched filter. Referring back to eq. (7-8), let us select 
S to be extremely small compared to the input noise samples. 
Under this condition, eq. (7-9) will become 

R nn .W = P (7-10) 

? 

where the non- zero term of Rgg:=(S ) is very small compared 
to E[x s x,.]. Recalling from matched filter derivation that 
the matched filter solution is given by 

R nn -W = C-S (7-11) 



where C is a constant and S is the 
case of point targets, both P and S 



signal vector, 
will be in the 



In the 
same form. 



P = [0 , . . . s 2 , . . . 0] 
S = [0 , . . . s , . . . 0] 



Since we can select the constant C arbitrarily, selecting 
C = S will cause eq. (7-11) to be identical to eq. (7-10). 

A practical implementation method to adapt the 
corresponding matched filter was performed in this work as 
the following: 

1. Select dj e S s 0 to satisfy equation (7-10). 

2. Since the constant C may be selected arbitrar 
ily and likewise, one of the filter coefficients can be 
selected arbitrarily, the estimation point coefficient was 
selected to be 1.00. Resulting simulated adaptation plots 
will be presented in this chapter. This modified Wiener 
filter coefficients adaptation curve, when compared to the 
theoretical matched filter coefficients, shows excellent agree 
ment . 
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To conclude the matched filter discussion, let 
us note that for the case of point target enhancement a 
modified adaptive Wiener filter will yield a corresponding 
matched filter. For targets other than point targets, in 
many practical cases we can find a corresponding matched 
filter counterpart. This issue is beyond the scope of this 
work and will not be pursued any further. 

e. The Selection of the Adaptive Loop Gain, y 

Appendix C, section C. presented few factors 
related to the selection of the loop gain. The first is the 
misadjus tment , given in eq. (C-63), and the second is the re- 
lated convergence time constant given in eq. . (C-64) . A more 
crucial parameter related to y, is the stability of the 
adaptive loop. 

In Chapter VIII we subjected the processing 
algorithm to the assumption that no a priori knowledge is 
available. Based on this assumption, an appropriate auto- 
matic gain control ("AGC") was needed. Various gain control 
loops were tested and the selected configuration is presented 
in Fig. 7-9. The "AGC" loop major design requirements were: 

1. Under all possible reasonable initialization 

of filter coefficients (W ) loop stability has to be assured 

[selection of W =0 should be avoided in any case) . 
o 

2. For any external selection of loop gain 
constant, C, the loop should not diverge (avoid overflow). 

3. The steady state misadj us tment , independent 
of background noise statistics, should be the same and should 
be controlled externally. 
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4. Subject to the above stated requirement, it 
is required to control the gain so as to yield the shortest 
adaptation time constant. 

The following measures were taken to satisfy the 
listed requirements: 

1. In order to avoid divergence of the filter 
output under any circumstances the output y(m,n) is limited 
to a selected maximal value. 

2. In order to assure a constant "misadj ustment" 
the total input power in each input sample set, has to be 
calculated to get the estimate of ( z — \ - p — ). This estimate 

llaCc 

is required since the misadjustment is given by: 



While estimating 
additional constr 



M = y • [trace* . 

(trace R xx ^ by A 1 * A 2 
aint is imposed on the 



a 

A 3 in Fig. 7- 9, 
trace measurement 



an 



G 2 ( k ) = i=i x i + I 



(7-12) 



As the filter adapts its coefficients, the output y becomes 

7 

very small compared to ^-]_ x i an< ^> after a limited number 

of adaption steps, if indeed convergence is achieved, 



G 2 (k) = 



£ Y • 

i= l X l 



As indicated 
externally . 
selected, by 



in Fig. 7-'9,- the gain constant C can be set 
Therefore, in extreme cases when a high C is 
mistake or intentionally, the second term in 
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eq. (7-12), namely,. |y(k)|, will assure the gain convergence, 

to prevent loop divergence. 

As the filter converges, the trace estimate 

becomes less noisy due to the smoothing operation done by . 

Eventually, after a few hundred adaption steps, the loop gain 

is stable and can be locked to its last value, as demonstrated 

in Fig. 7-9 by the memory, A^, and the switch, SW1. 

By utilizing the "AGC" configuration presented 

in Fig. 7- 9, the misadjustment is indeed constant under all 

circumstances since the loop gain p = C - — ^ — and the mis- 

tr k xx 

adjustment given by ( C - 63) is M = ytr R xx = C (constant). 

Then, in Appendix C, section C, the adaption time constant 
was calculated as well and was given by (C-67) (for equal 
eigenvalues) . 

I = J_ 

T mse 4y tr R xx 4C 

For unequal eigenvalues, x is also related to i but in 
a more complex form. 

Therefore, this "AGC” configuration could yield 
the required performances, for a proper selection of the 
constant C. For example, a misadjustment of 10% can be 
assured by selection of C = yjy . Fixing C will result in a 
given adaptation time. 

In order to accelerate the adaptation period, a 
gain switching option was tested such as demonstrated in Fig. 

7- 9 ()b) . In the gain switching case, an initial gain constant, 
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as high as possible, is selected to achieve a fast conver- 
gence toward steady state. After a certain number of adap- 
tation steps, the gain constant is switched to a low value to 
assure the proper required misadj us tment . When some a priori 
knowledge is available about background noise statistics, the 
gain switching point can further optimize the performance. 

There is a clear trade-off that limits the util- 
ization of the gain switching method. Whenever the initial 
guess, W q (the initial filter coefficient) is close to the 
steady state solution, applying the gain switching technique 
will not be fruitful and sometimes will even be the worst. 

On the other hand, when initial guess of W is assumed to be 
far from the steady state solution, the gain switching tech- 
nique will indeed shorten the adaptation time. 

f. The Required Number of Independent Pixels to 

to Achieve Adaptation 

Throughout all the preceding discussion the image 
statistics were assumed stationary. If one assumes that the 
stationarity is confined to limited image subzones, the 
adaptive filter input samples should be confined accordingly 
to a subzone. The question raised is, what is the minimum 
number of independent pixels sufficient to achieve filter 
adaptation? Intensive computer simulation indicated that 
for a 3x5 filter size a subimage of size 6x6 is sufficient 
to adapt filter coefficients. It was noticed that adaptation 
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rate is slowed down as the number of independent pixels i's 
reduced. Simulation results for small size image filter 
adaptation will be presented at the end of this chapter. 



E. COMPUTER SIMULATION RESULTS 

1 . Introduction 

In this section typical computer simulation results 
will be presented. The type of filter selected to be presented 
is a spatial 3x3 size filter, optimized to enhance point 
targets and suppress correlated noise with variance = 4.0. 
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The following plots will be shown: 



1. The filter coefficients adaptation plots, compared to 
the theoretical coefficients. 

2. The learning curve (MSE as a function of adaptation 
steps) . 

3. The gradient estimate curve as a function of adaptation 
steps . 

The simulation procedure was conducted in the following way: 

1. Generate the correlated noise image with specified 

statistical properties. 

2. Adapt the filter coefficients. 

3. Design a theoretical Wiener or matched filter based on 
the known statistical properties. 

4. Plot the coefficients adaptation curve compared to the 
theoretically designed coefficients. 

Three sets of adaption plots will be presented: 

1. Adaption of the spatial filter to an image with first 
order Markov model statistics. 

2. Adaption curves for second order Markov model clutter. 

3. Small size image adaptation curves. 
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2 . Adaption to First Order Markov Model 



Fig. 7-10 presents the unprocessed clutter image 
generated by first order Markov process with correlation 
factor 

P H = P V = 0,99 

The theoretical filter coefficient and the steady state 
adaptive coefficients are shown as well in Fig. 7-10. 

Fig. 7- 11 (a) , (b) , (c) presents the filter coefficients 
adaptation curves as a function of adaptation steps. The 
theoretical coefficients level is presented as well in each 
adaptation curve and it is seen that the adaptive filters 
coefficients come indeed very close to the theoretical value. 

Fig. 7-12 is the learning curve and it shows the MSE 
at filter output as a function of adaptation steps. 

Fig. 7-15 shows the sum of the absolute values of all 
coefficients correction signals. 

3 . Adaptation to Second Order Markov Model 

Fig. 7-11 presents the unprocessed clutter image 
generated by second order Markov process with correlation 
factor 

Pj_j — Py ~ 0.95 

and spatial frequency 

B H = (3 y = 0.628 

The theoretical and the steady state adaptive filter coeffi- 
cients are shown as well in Fig. 7-14. 

Fig. 7 - 1 5 (a) , (b) , (c) are the filter coefficients 
adaptation curves. Fig. 7-16 is the learning curve and 
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The Unprocessed Image 
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Figure 7-10 Unprocessed Image and Filter Orientation 
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Figure 7-14 Unprocessed Image and Filter Orientation 
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Fig. 7-17 is the sum of the absolute values of all coeffi- 
cients correction signals. 

4 . Adaptation to Small Size Image 

Fig. 7-13 is the unprocessed first order Markov model 
clutter image. A section of 6x6 pixels is selected as shown 
in Fig. 7-18 to adapt the filter coefficients. The adaptive 
process is accomplished within this small site section. The 
resulting filter coefficients adaptation curves are shown in 
Fig. 7-19 (a) , (b) , (c) . The learning curve is shown in Fig. 
7-20 and the sum of coefficients correction signals is shown 
in Fig . 7 - 21 . 

5 . Simulation Results Conclusions 

The major computer simulation results conclusions are 
described below. 

- The adaptation rate decreases as correlation 
coefficient increases. 

- An impirical relation was approximated to indicate 
the total number of adaptation steps N1 (total) , as a function 
of p H or p v and S H or gy. 

- For the first order Markov process: 

NX (total) « k x [ 1 . 0 - (.P max ) 2 ] 
where p indicates the higher correlation coefficient 
(between p^ and Py) • Typically, k^ = 250 - 500, depending on 
loop gain constant. 

- For second order Markov process, the following 
approximated empirical formula can be used: 
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N2Ctotal) - N1 (total) *[ (1 + 1/cos 8 )] 

m 3.X 

8 is the higher 8. 

- In general, it was noticed that the adaptation 
rate is slower in the second order Markov process than in 
the first order Markov process. 

- The adaption time of a small sized image is slower 
compared to larged sized image (by large image it is meant 
that the image size is very large compared to the search box 
size) . 

- A complementary study with regard to suboptimal 
filters was carried out in the study phase presented in 
Chapter VI. One of the conclusions there was that a sub- 
optimal filter designed to operate in the presence of highly 
correlated clutter is the preferred suboptimal filter. Com- 
bining this conclusion with the fact that adaptation time is 
slow in case of high correlation - it is always recommended 
to select the filter initial value as if it were designed to 
operate with the highest correlation factor . 

F . SUMMARY 

In this chapter a spatial adaptive filter concept was 
discussed and implemented. The adaptation method presented 
in the chapter is applicable for enhancement of targets of a 
given shape cluttered hy highly correlated noise. 

In addition, a practical method to modify a special Wiener 
filter - intended to detect point targets - into a matched 
filter was presented. This method enables a direct adaptation 
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of a matched filter in the case of point target. The 
suggested method was tested against simulated images and 
demonstrated the expected performance - namely, the adapta- 
tion of the filters into a steady state optimal solution. 

Various practical implementation issues were mentioned 
in the chapter such as: 

- The image scanning pattern 

- The "desired response" derivation 

- The adaptive loop gain setting 

- The typical adaptation rates 

- The required number of independent pixels. 

This last issue will be utilized in the next chapter 
where nonstationary statistics will be assumed. 
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VIII. ADAPTIVE SPATIAL FILTER II - 

NONSTATIONARY CLUTTER STATISTICS 



A. INTRODUCTION 

In this chapter, an adaptive spatial filter technique 
capable of processing images having nonstationary clutter 
will be developed. Its basic idea is described in the fol- 
lowing. More detailed explanations will be given in following 
sections . 

The first assumption is that an image of nonstationary 
clutter statistics can be divided into subimages within which 
the clutter statistics are stationary. 

A spatial filter will be assigned to each of these sub- 
images. An adaptive algorithm is then applied to adapt the 
assigned filter to its optimal performance. An automatic 
procedure is applied to terminate the adaptation. The optimal 
filter is then applied to the subimage. 

Following this, an optimal threshold will be determined 
and applied to the adaptively filtered subimage for target 
detection. 

B. THEORETICAL DEVELOPMENT AND PRACTICAL CONSIDERATIONS 

1 . Subimage Division 

In this thesis, a procedure of dividing the image 
into square shaped small subimages is developed. From the 
processing gain requirements, the smallest subimage size is 
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preferred hecause the chance is better to have "stationary" 
statistics in a small area. This makes possible the design 
of an optimal spatial filter. On the other hand, smaller 
subimage size requires more adaptive filters to cover a given 
image size and increases the computation load. A trade-off 
should be made for an optimum subimage size. The smallest 
size used is 6x6. Using this small size and simulated images 
containing Markov model correlated noise, it was found that 
the adaptive filter was able to converge to the optimal value 
in all cases tested. 

2 . Adaptive Spatial Filter 

For each subimage, adaptive spatial filters developed 
in the last chapter can be used because we assume that the 
statistical properties in the subimage are stationary. There- 
fore, detailed explanation for the adaptive spatial filter 
will not be repeated here. Instead, only the important steps 
will be highlighted in the following. 

a. Initialization of the Filter Coefficients 

It can be started in two ways. We can assume 
a statistical characteristic and design a spatial filter 
accordingly. Alternatively, we can calculate the statis- 
tical characteristics of the image by ignoring the nonsta- 
tionarity which might exist in the image and design a spatial 
filter accordingly and use it as the initial filter. 
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b. Selection of Adaptive Loop Gain 

It was discussed in detail in Chapter VII, section 
G.4. In our computer simulation, a loop gain of 0.4 was used 
most frequently. 

c. Termination of Adaptation 

The procedure to terminate the adaptation was 
discussed in Chapter VII, section G.2. In our computer simu- 
lation, the following parameters are used: 

n = 36, N = 360, F = 5 % 

d. Adaptation Time 

It depends on the following factors: 

Filter initialization 

Adaptive loop gain 

Statistical characteristics of the image. 

The last factor is beyond the control of the filter designer. 
However, design of the initialization and loop gain could have 
a significant effect on the adaptation time. It was found, 
from computer simulation, that using a loop gain in the range 
of j to it took at most several thousand adaptive steps to 
reach the optimal filter in most cases. 

C. COMPUTER SIMULATION RESULTS 

1 . Generation of Nonstationary Images 

The nonstationary image used in this chapter is the 
same one introduced in Chapter VII, section G.4. It is an 
image of size 64x64 pixels, and consists of four 
quadrants, each of which has its own statistical properties. 
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They are separately shown in Figures 8-1 through 8-4, 
which include their statistical properties and optimal spatial 
filters together with the theoretical processing gain of these 
filters . 



Fig. 8-1 Upper left quadrant 
Fig. 8-2 Upper right quadrant 
Fig. 8-3 Lower left quadrant 
Fig. 8-4 Lower right quadrant 

The composite nonstationary image is shown in Fig. 8-5(a) 
by a 9 gray level printout and in 8-5 C b ) by a 3-D plot. 

2 . Image Subdivision and Filter Adaptation 

The image is divided into square shaped subimages of 
size 8x8 pixels, as shown in Fig. 8-6(a). An adaptive filter 
was assigned to each subimage and was adaptively changed until 
it reached the optimal performance. In Fig. 8-6 (b) the adaptive 
filter coefficients of several selected subimages are shown and 
can be compared with the theoretical optimal filter. The 
resulting processing gain and optimal threshold are shown in 8-6(c). 
The threshold was determined by assuming Gaussian probability 
density function, and for a probability of false alarm in the 
range from 1 x 10 _1 <_ Pp^ < 5x 10 J 

3 . Comparison of Filter Performance 

The following figures of merits are used to evaluate 
the filter performance. 

For clutter suppression: 

a. Processing gain (P.G.) 
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Figure H-u(a) Computer Print and Statistical Properties oi Upper Left Quadrant Image 

and its Optimally Designed Spatial Matched Filters. 



3-D Plot of Upper Left Quadrant 
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' igure 8-2 (a) Computer Print and Statistical Properties of Upper Right Quadrant 
Image and Its Optimally Designed Spatial Matched Filters. 




Figure 8-2(b) Three Dimensional Plot of the Computer 

Generated Upper Right Quadrant Image. 
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iguie 8 3(a) Computer Print and Statistical Properties of Cower Left Quadrant 
Image and Its Optimally designed Spatial Matched Filters. 



3-D Plot of Lower Left Quadrant 



Z 





Figure 8-5(b) 



Three Dimensional Plot of Computer 
Generated Lower Left Quadrant Image. 
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Figure 8-4(a) Computer Print and Statistical Properties of Lower Right Quadrant 

Image and Its Optimally Designed Spatial Matched Filters. 




Figure 8-4(b) 



Three Dimensional Plot of Computer 
Generated Lower Right Quadrant Image 
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3rd 



4th 



Figure 3-5(a) Nine Gray-level Computer Print of Four Quadrant 
"Nonstationary” Simulated Clutter Image, 







Figure 8-5(b) Three Dimensional Plot of "Four Quadrant" 
"Nonstationary" Simulated Clutter Image. 
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The Image Division to Many 
Subimage Zones 



Image Image Image Image Image 

Zone 1 Zone 2 Zone 3 Zone 4 Zone J 

FILTER I FILTER 2 FILTER 3 FILTER 4 FILTER J 




Figure 8-6(a) Division of a Given Image Into Subimage 
Zones for Performing Filter Adaptation. 
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3rd Quadrant 



4th Quadrant 



Figure 8-6(b) Adapted Spatial Filter Coefficients for Each 

8 by 8 Pixels Adaptation Subimage Zone. 
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Processing Gain and Threshold Level in Each Subimage 




Estimation point in the inner sectors is in the 
search box center. 

Figure 8-6(c) Processing Gain and Threshold Level for 

Each 8 by 8 Pixels Adaptation Sub image Zone. 
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For target detection: 

b. Probability of detection (P^) 

c. Probability of false alarm (P CA ) 

Sets of comparisons will be presented, 
a. For Background Clutter Suppression 

The first set is to compare adaptive and non- 
adaptive filters in terms of background clutter suppression 
and probability of false alarm. It consists of five figures. 

Fig. 8-7 shows the adaptively filtered and thresh- 
olded background clutter image, using a subimage size of 
8x8 pixels. The resulting probability of false alarm is also 
listed in Fig. 8-7. 

Fig. 8-8 shows nonadaptively filtered and thresh- 
olded background clutter image and consists of four figures. 

Fig. 8-8(a) is the result using the optimal spatial filter 
designed for the first quadrant. Figures 8-8 (b), (c) and (d) 

are the corresponding results using spatial filters optimally 
designed for the second, third, and fourth quadrants respectively. 
Probabilities of false alarm for each case are also listed in 
these figures. 

Comparing Fig. 8-7 and Fig. 8-8 revealed that the 
probability of false alarm of the adaptive filter case is 
about '22 times lower than that of the nonadaptive filter case. 

The second set of comparisons is similar to the 
first set for comparing clutter suppression and Pp^. However, 
the adaptive spatial filter in this set used a large subimage 
size of 16x16 pixels. Only three figures are used. 
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Figure 8-7 Computer Print of False Alarms of Adaptively 
Filtered and Thresholded Four Quadrant 
"Nonstationary" Simulated Clutter Image. 
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PROBABILITY CF FALSE ALARM IS 0.1431 



Figure 8-8 (a) Computer Print of False Alarms of Nonadapt ively 

Filtered and Thresholded Four Quadrant "Non- 
stationary" Simulated Clutter Image Using Spatial 
Filter Optimally Designed for First Quadrant. 
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PROBABILITY OF FALSE ALARM IS 0.1462 



Figure 8-8 (b) Computer Print of False Alarms of Nonadaptively 

Filtered and Thresholded Four Quadrant "Non- 
stationary" Simulated Clutter Image Using Spatial 
Filter Optimally Designed for Second Quadrant. 
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PROBABILITY OF FALSE ALARM IS 0.1596 



Figure 8-8(Y:) Computer Print of False Alarms of Nonadaptively 

Filtered and Thresholded Four Quadrant "Non- 
stationary" Simulated Clutter Image Using Spatial 
Filter Optimally Designed for Third Quadrant. 
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PROBABILITY OF FALSE ALARM 


IS 


0.1335 











Figure 8-8(d) Computer Print of False Alarms of Nonadaptively 
Filtered and Thresholded Four Quadrant "Non- 
stationary" Simulated Clutter Image Using Spatial 
Filter Designed Optimally for Fourth Quadrant. 
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Figure 8-9 shows the adaptively filtered and 
thresholded clutter image and its P~ , . 

Figure 8-10 compares the performance of adaptive 
filters using 8x8 and 16x16 subimage areas. It can be seen 
that a smaller subimage size offers better results. Figure 
8-10 indicates that in 7 out of 16 (431) subimage areas of 
size 16x16 a transition in the statistical properties occurs 
inside the subimage area. On the other hand, Fig. 8-6(c) 
indicated that in the 8x8 subimage areas the transition in 
the statistical properties occurs only in 15 out of 64 (25%). 
This difference explains the improvement in the Pp^ as sub- 
image areas size is reduced. 

Figure 8-11 shows the nonadaptive filtered and 
thresholded image using a spatial filter optimally designed 
for the lower right quadrant. The resulting Pp^ is also 
listed in the figure. 

Comparing figures 8-9 and 8-11 revealed that the 
Pp^ of the adaptive filter case is about five times lower 
than that of the nonadaptive filter case. The difference is 
not as big as a factor of 22 for the smaller subimage size 
case shown in figures 8-7 and 8-8. 

For Clutter Suppression and Target Detection 
The third set of comparisons is for both clutter 
etection. For this purpose, an en- 
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THE ADAPTIVELY FILTERED IMAGE 
PROBABILITY OF FALSE ALARM IS 0.0078 
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Figure 8-9 Computer Print of False Alarms of Adaptively 
Filtered and Thresholded Four Quadrant 
"Nonstationary" Simulated Clutter Image. 
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Figure 8-10 Processing Gain and Threshold Level for Each 
16 by 16 Pixels Adaptation Subimage Zone. 
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PROBABILITY OF FALSE ALARM IS 0.0383 

Figure 8-11 Computer Print of False Alarms of Nonadap tive ly 
Filtered and Thresholded Four Quadrant "Non- 
stationary" Simulated Clutter Image Using Spatial 
Filter Designed Optimally for Fourth Quadrant. 
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THE ORIGINAL TARGETS LK4TI0N 

TARGETS LEVEL IS 1.50000 TIMES THRESHOLD LEVEL 



Figure 8-12 Locations of Point Targets Added to the Four 
Quadrant "Nonstationary" Simulated Clutter 
Image for Calculating Probability of Detection. 
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TARGETS TOGETHER WITH BACKGROUND NOISE 






Figure 8-13 Computer Print of Four Quadrant "Nonstationary" 



Simulated Clutter 



Image 



Including Added Point 



Targets 



Shown in Figure 8-12 



•ns 



The target intensity is proportional to the threshold level 
at each subimage zone. 

Figure 8- 14 presents the adaptively filtered and 
thresholded image in which all false alarms have been removed 
so that the probability of detection can be calculated. It 
is included in the figure. 

Figure 8-15 (a) , (b) , (c) , (d) present the nonadaptively 
filtered and thresholded image using spatial filters designed 
optimally for first, second, third and fourth quadrants respec- 
tively. Resulting probability of detection for each case is 
included in corresponding figures. 
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THE ADAPT IV ELY FILTERED IMAGE 

Figure 8-14 Computer Print of Adaptively Filtered and 
Thresholded Four Quadrant "Nonstationary" 
Simulated Clutter Image Including Added 
Point Targets. 
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PROBABILITY OF DETECTION IS 0.731 

Figure 8-15(a) Computer Print of Nonadaptively Filtered and 

Thresnolded Four Quadrant "Nonstationary" 
Simulated Clutter Image. Spatial Filter Used 
was Designed Optimally for the First Quadrant. 
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Figure 8- 15(b) Computer Print of Nonadaptively Filtered and 

Thresholded Four Quadrant "Nonstationary" 
Simulated Clutter Image. Spatial Filter Used 
was Designed Optimally for the Second Quadrant. 
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PROBABILITY OF DETECTION IS 0.694 

Figure 8-15(c) Computer Print of Nonadapt ively Filtered and 

Thresholded Four Quadrant "Nonstationary" 
Simulated Clutter Image. Spatial Filter Used 
was Designed Optimally for the Third Quadrant 
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Figure 8- 15(d) Thresholded Four Quadrant "Nonstationary" 

Simulated Clutter Image. Spatial Filter Used 
was Designed Optimally for the Fourth Quadrant. 
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D . SUMMARY 



An adaptive spatial filter algorithm based on least mean 
square criteria and using a steepest descent technique is 
shown to be effective in suppressing nonstationary clutter 
based on computer simulation of Markov model clutter noise. 

The probability of false alarm of the adaptive filter 
case is lower than that of nonadaptive filter case by 5 to 
22 times for the examples used in this study. 

The probability of detection of the adaptive case is 
slightly better than that of the nonadaptive case when the 
threshold level is kept equal for both cases. It must be 
emphasized that Pp^ for the nonadaptive case is 5 to 22 times 
higher. On the other hand, if a CFAR detection is applied 
by adjusting the threshold for the nonadaptive case, severe 
reduction of P, will occur. 
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IX. SUMMARY AND CONCLUSIONS 



A. SUMMARY 

1 . Introduction 

New infrared surveillance and guidance systems using 
large mosaic detector array focal planes require the detec- 
tion of targets deeply buried in background clutter noise by 
several tens of db ' s . New image processing techniques must 
be developed to accomplish these heretofore unachievable goals. 
They must also be relatively straightforward so that on-board 
implementation using moderate LSI signal processors is pos- 
sible to meet real time processing requirements. They should 
also be flexible enough so that a focal plane processing 
algorithm is able to adapt itself to changing background scenes. 

This thesis develops five focal plane processing pro- 
cedures which consist of nonrecursive statistical filters for 
background clutter suppression followed by thresholding to 
initiate the target detection. Nonrecursive approaches are 
taken because they make extensive use of both the spatial and 
temporal behavior of the image signal in the neighborhood of 
the pixel of interest. On one frame of the image, signals 
of pixels all around the pixel of interest are used. In 
addition, signals at the same pixel location from other frames 
before and after the frame of interest are also used in the 
focal plane processing. New integrated circuit tapped delay 
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line devices may be used effectively for their hardware 
implementation. Statistical approaches are used because it 
is believed that a. collective statistical description can be 
found for different background scenes if they belong to the 
same class of terrain, weather conditions, etc. Therefore, 
it is possible to design a statistical filter which is able 
to estimate a weak target buried in noise and is not limited 
to only a few background scenes. 

The concept of using nonrecursive statistical filters 
for background clutter suppression was originated by Lt. D. 
Bar-Yehoshua [14] in his development of nonadaptive spatial 
filters based on the minimisation of mean square error. This 
thesis made extensive extension and more thorough developments 
of four nonadaptive filters in both spatial and temporal do- 
mains. In addition, a new adaptive spatial filter is developed. 

2 . Nonadaptive Focal Plane Processing 

In Chapter II, the basic concept of the nonrecursive 
statistical filter is expanded and generalized to include all 
nonadaptive and adaptive, spatial, temporal and mixed spatial- 
temporal or temporal -spatial filters. A unified processing 
concept is presented for clutter suppression based on the 
statistical spatial and temporal behaviors of the image in 
the neighborhood of the pixel of interest. Two statistical 
criteria for filter design are introduced. One is based on 
the minimization of mean square error criterion which leads 
to the MMSE filter or Wiener filter. The other is based on 
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the maximization of signal to noise ratio criterion which 
leads to the statistical matched filter. 

Detailed developments of the principles, configura- 
tions, and design procedures for the filters are developed 
in Chapter III for spatial filters using a single image 
frame and in Chapter IV for sequential spat ial - temporal 
filters using several frames of images. These three dimen- 
sional statistical filters are useful in situations when 
multiple frames of images are available, either when the 
surveillance and guidance systems are already in their track- 
ing modes after the targets have been acquired or when the 
fields of view of the systems are stabilized to cover an 
area for a number of frames . 

In Chapter V, the thresholding procedure for a two 
dimensional image is developed for the initialization of 
target detection. Two CFAR (constant false alarm rate) 
procedures to set the threshold level are developed. One 
is for Gaussian probability density functions. The other 
is for the more general non-Gaussian PDF. 

In Chapter VI, the background clutter suppression 
filters and thresholding are combined together in four non- 
adaptive focal plane processing procedures. Their properties 
and performance are studied in detail, using two types of 
images. One type is a set of twenty frames of real world 
infrared images with 10% uniform drift between frames. The 
second type is computer generated images which use two Markov 
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processes to model the correlated noise which is assumed to 
be a good simulation of the real world background clutter. 

It has recently been brought to our attention that the Markov 
model was found to be an acceptable simulation of the back- 
ground clutter based on recent work on a tactical infrared 
hostile weapon location system.* With the computer generated 
model, a wide variety of images with varying statistical 
properties can be generated with ease. 

Both background clutter suppression performance and 
target detection performance are evaluated by computer sim- 
ulation using both types of images. Using the real world 
infrared images, a background clutter suppression of 27 db 
is achieved by a single frame spatial filter. Using a five 
frame spatial- temporal filter, a suppression of 87 db is 
achieved. These results are the most important testimony 
of the practical usefulness of the nonadaptive statistical 
filters because they have been used by other research groups 
and considered as one of the "standard" test scenes. Com- 
parison with results of other groups is difficult in this 
unclassified study because the absolute calibration of the 
image signal is a classified number and was not made avail- 
able. It should be acknowledged that the real world infrared 
image test described by figures 6-11 to 6-18 was made by 
Captain D. Hilmers and was reported originally in Ref. 50. 

*J. Otazo, MIT Lincoln Laboratory, private communication, 
and to be presented in the SPIE 1979 symposium, Washington, 
D.C. , April 17-18, 1979. 



237 



3 . Adaptive Focal Plane Processing 



In Chapters VII and VIII, adaptive statistical spatial 
filters based on the minimization of mean square error cri- 
teria are developed for clutter suppression. In Chapter VII, 
the adaptive procedure is deyeloped using computer generated 
"stationary" images whose statistical spatial correlation 
properties are uniform. In the following chapter, the pro- 
cedure is applied to computer generated "nonstationary” images 
which have four different statistical spatial correlation 
properties in a given image. It is found that using a small 
adaptive array of 8x8 pixels, better clutter suppression can 
be achieved compared with nonadaptive statistical spatial 
filters. Using the adaptive approach, the resulting proba- 
bility of false alarm (Pp^) of the adaptively filtered and 
thresholded image was found to be more than twenty times 
better than that of the four possible nonadaptive filters. 

B. CONCLUSION 

In conclusion, five nonrecursive statistical focal plane 
processors have been developed. The usefulness of this new 
concept of clutter suppression has been demonstrated by apply- 
ing these nonadaptive filters to a set of real world infrared 
images. Although the adaptive spatial filters have been 
tested on computer generated images, it is believed that good 
clutter suppression can also be achieved on real world images. 
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Neither the concepts of nonrecursive filters nor that of 
statistical estimation are new. However, the combination of 
these two concepts to process two and three dimensional images 
of very low signal to noise ratio for suppression of "colored" 
clutter noise and for substantial enhancement of target signal 
to noise ratio up to several ten's of db ' s has not been 
reported before. The contributions of this thesis to the 
advancement of focal plane processing techniques for the 
enhancement of point/line targets and the suppression of 
clutter are highlighted in the following. 

1 . Nonadaptive Focal Plane Processing 

a. One dimensional statistical estimation techniques 
based on the minimization of mean square error criterion 
(Wiener filter) and the maximization of signal to noise ratio 
criterion (matched filter) are extended to two dimensions for 
a single frame and to three dimensions for multiple frame 
focal plane processing. They are different from other mul- 
tiple dimensional image processing studies because special 
emphasis is given to the suppression of clutter which is 
modeled by correlated noise. Furthermore, target signals 
several ten's of db ' s helow clutter noise are considered. 
Images of such low signal to noise ratio have not been con- 
sidered until motivated by the recent interest in smart sensor 
systems . 

b. The clutter suppression is accomplished by per- 
forming an estimate of the target. This is different from 
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other approaches which first estimate the clutter and then 
suppress it by subtraction from the estimate of the target signal. 

2 . Adaptive Focal Plane Processing 

a. One dimensional nonrecursive adaptive filter 
techniques based on the least mean square criterion are ex- 
tended to two dimensional single frame focal plane processing 
with special emphasis on the suppression of correlated clutter 
noise several ten's of db's higher than the signal of interest. 

b. A procedure to derive the desired response corre- 
sponding to different types of target signals is developed 
which is used in the adaptive process as the optimal response. 

c. The adaptive loop gain is set by an "automatic 
gain control" procedure which includes several requirements 
simultaneously to yield constant misadjustment for varying 
statistical properties of the clutter, to prevent overflow 
and to assure adaptive convergence. 

d. A procedure for dividing images whose statistical 
properties are nonstationary and unknown a priori, into sub- 
image zones for filter adaptation is developed which minimizes 
the degree of nons tationarity and improves the performance of 
adap tation . 



240 



APPENDIX A [Ref. 14] 



GENERATION (MODELING) OF TWO-DIMENSIONAL RANDOM FIELDS* 

A. INTRODUCTION 

Any monochromatic image can be modeled by specifying its 
value (gray level) f(m,n) at each spatial coordinate (m,n). 

An ensemble of such images can be modeled by interpreting 
f(m,n) as a random field. A class of images with particular 
spatial correlation can be modeled by a random field with 
the same autocorrelation function. 

When spectral factorization in two dimensions is possible 
[Ref. 56], the autocorrelation function of a random field is 
related to the transfer function of a system that generates 
the random field from a set of uncorrelated variates. Then, 
specifying a class of images by their autocorrelation func- 
tion specifies a dynamic model representation for these images. 

For a discrete random- f ield , the dynamic model is a 
partial difference equation that can be represented by a two 
dimensional recursive relation. The problem in this Appendix 
is to find a linear equation that represents a given random 
field. A discrete two dimensional random process with two 
kinds of autocorrelation functions will be considered. 

The presented method uses the spatial frequency domain 
to generate the requested random field. 

*Excerpt from thesis by D. Bar-Yehoshua, "Two Dimensional 
Nonrecursive Filter for Estimation and Detection of Targets," 
Naval Postgraduate School, June 1974 [14]. 
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B. SPATIAL FREQUENCY DOMAIN METHOD (TRANSFER FUNCTION) 

1 . Introduction 

The discrete power transfer relation, for linear 
filtering of random signals, is useful for designing signal 
processors. If we like to generate a random field with given 
autocorrelation function, we can find a transfer function 
through which an uncorrelated random field at the input (white 
noise), will give a defined discrete power spectrum density at 
the output. 

Figure A-l shows the relationship between a white 
noise input signal passing through a linear filter and the 
resulting output. In this case, the discrete power density 
of the input is the 2-D unity function, and the discrete 
power density of the output equal to the absolute square of 
the filter transfer function. 

The transfer function method procedure can be written 
as follows: 

a. From the autocorrelation function of the random 
field to be modeled, find the discrete power spectral density 
using z transform. 

s Ul,- 2 ) = 2{R(m,n) } (A-l) 

b. Assume that the discrete power spectral density 

SCz x z 2 ) = |H(z 1 > z 2 )| 2 = H(z 1 z 2 )*H(z 1 ' 1 Z2 -1 ) (A-2) 



2 42 



is : 



Deterministic 

Input 



04 

M 



CNJ 



fv| 

v — / 

tu 



04 

M 



t4 

X 

II 




c 

£: 

' — / 

m 



4 -« 

CO 



04 



04 

M 



M 

X 



04 



(A 

C 

<D 

H 3 



03 

Jr 

4 -> 

U 

<D 

X 

(/) 

Jr 

<D 

£ 

O 

a* 




c 

•K. 

s 

V_/ 

X 



04 

M 



M 

X 

II 




4-1 



04 





/ \ 

04 


/ N 

a 


Nl 


<L> 

+-> 


11 


/ N 

04 




/ N 


M 


•s. 


*— 1 


' 4> 


/ — n 


M 




c 




6 


tsl 


Jr 


c 




c 




rH 


v — / 


v — ' 


U 




rH 




6 


t4 


DO 


DO 


to 


JS 


K ^ 


s 


v — ' / 




DO 4 


\ DO 


♦ rH 






V J l 


DO 


o 


X 


CO 


X 

' / 


00 


o 


DO 



04 

N 



M 

X 



04 

f 4 






CO 



DO 

OO 







t 


7 


/ — \ 






rH 


/ — \ 


/ N 


>-* 


c 


0 ) c 




** 


to 


s 


e 


rH **H £ 


v — ' 


V — / 


o 



03 



4h 


(J 


II 


/ s 


4> II 


/ \ 


X 






04 


4-» 


O 




4-> 


/ — \ 


t4 


•rH / — \ 


M 




to 


c 




X c 






•H 




rH 


5 - 


r*“ 




c 




t4 


£ 


M 








S J 




^ — / 




*rH 


4h 


X 


_ 4h 


CO 


s 


6 






E 




O 4J 


u 


+-» 




O 4-> 




03 3 


<D 


3 




T3 3 




C X 


4~> 


X 




^ c X 




03 C 


0) 


C 




X 03 C 




X nr 


X 


nr 




v — ■ / X nr 





243 



Figure A-l. A filter with deterministic and random inputs 
(a) a general case, (b) with impulse and white 
noise inputs. 



Try to separate S(z^ z 7 ) in order to get H(z^ z 9 ) as a trans- 
fer function of the random field model, according to eq. A- 2 . 

c. Using an uncorrelated random field as an input to 
H(z^ Zj ) , find the dynamic model of the designed random field 
in the z domain. 



GCZi z 2 ) = H C Z 1} z 2 ) * u Cz 1 z 2 ) (A- 3) 

where u(z, z ? ) is a random white noise signal. 

d. Convert the dynamic model to a spatial model by 
inverse z transform. In this process it is assumed that the 
random field is stationary and ergodic. In this section the 
dynamic model for two autocorrelation functions will be con- 
s idered . 

2. Dynamic Model of One-Dimensional 7 ( . 

Autocorrelation Function R(t) = c^ r p 1 T 

Given the autocorrelation function of one-dimensional 
random signal as: 

R(m) = o R ?p |m! (A-4) 



The discrete power spectral density will be: 

g r 2 * (1-P 2 ) 

S(z) = — ^ 



(l-pz)*(l-pz ) 



(A - 5) 



J T 

°R‘ 1 ' p ' 



r, z 

Or- 1-0 



l-pz 1 - pz 

So, the transfer function will be 



-1 



H(z) H(z _1 ) 



H(z) = 



Or- 



/ 



1 - 



pz 



(A- 6) 
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and the output Lor such transfer function: 



gU) = H ' - ) • J' ( z ) 




1_ PZ ' 1 



F(z) 



If we replace the input by white noise, we will get 
r ae dynamic model in z domain: 



gU) = pz'-g(z) + a^fl-p 2 'u(z) ( A- 7 ) 

or in time domain: 

g(n) = pg(n-l) + a^fi-p -tiCn) (A-8 ) 



where u(nl is a white noise with variance = 1. 

Eq . (A-8 ) represents a dynamic model for a one- dimens iona 1 

first order Markov process. 

3 . Dynamic Model of One-dimensional Autocelleration 
Function R(n) = Gp &p 1 n '.cos ( <u» nj 

Given the autocorrelation function: 

R(n) = a R ?p ^ n [cos(wnj ( A-9 ) 



The discrete power spectral density will be: 

a + bz 

.R-./l - 



S(z) = a D 2 .(l-p 2 )*- 



2 2 

(1- 2p*cos to z +■ p“z ) 



a + b: 



- 1 



— 2 - 2 

(1- 2p*cos uJz + P ‘Z 



(A- 1 0) 



where 



a = 



b = 



fj-2 p'CCSai + p 2 ' + ]^1 + 2p*cos oi ~p~ 

2 



f l 'l-3z» c osai+p - fl + 2p cos u> + p 



(A- 11) 



Note: 1 - 2 c o s w * >0 for 0 < o < 1 . 
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( A- 1 2 ) 



Separating S(zJ to ll(z)*H(z 1 ) , we get: 



H(z) = o R ./ 1-p 2 ' • 



a + bz - 1 



1 - 2p'COS W Z + p ♦ z 



There are two ways to represent the dynamic model of 
such transfer function. 

a. First way: 

G(z) = H(z) • F(z) 

Substitute H(z) from eq . (A-12) and represent f as a white 

noise, we will get: 




The dynamic model in time domain will be: 




•( a* u (n) + b • u (n- 1 ) ] 

( A- 14) 



where u(n) is a random white noise. 



The initial condition will be: 



assume 



g(-i) * g C 03 = u ( o) = o 



we will get: 



g ( 1 ) * o R -u(l) 



g ( 2 ) = 2 p'Cosu • g ( 1 ) + o R ^l-p »[au ( 2) + bu ( 1) ] . 



(A - 1 5 ) 



b. Second way - state variable representation 



Let us define: 




CA-16) 



,4 6 



so that: 



G(z) 



H (. z ) * z 1 * F 1 ( z ) 



- 1 



U) 



1 - 2p*cos luz 



- 1 




(A - 1 7) 



G (z) = 2* P * cos u) • z •G(z)-p.“z'. 2 G(z)+z <1 'F'(z) (A -18) 
Let us define: 

p 1 (z) 2 G ( z) 

p 2 (z) £ P 2, z' X »G(z) (A - 1 9 ) 

From eq. (A -18) and eq. (A-19), we get: 

p x (z) = 2p*cos uj-z ^ P -^ ( z ) - z _1 .p 2 (z) + z _1 *f' (z) 
2 - 1 

p 2 (z) = P. z • P 1 ( z ) 

Its expression in the time domain is: 

p^(n) = 2p -cosw*p^(n- 1) - p ? (n-l) + f '(n-l) 

2 

p 2 (n) = P • p 1 (n-l) (A-20) 



Expressing eq. (A-20) in the state variable form, we have: 



|~ P 1 (n) 




2p»cos w i - 1 
7 * 


' 

* 


«• — 

P x (n- 1) 


f 

1 

►— 1 

J 


'i_P2 (n ^ 




P 4 i 0 

. 1 




P 2 (n-l) 


1 

O 



for n > 1 



where the output is defined as: 



(n-l) 



(A- 21) 



g(n) = 



0] • 



P x (n) 
P 2 (n) 



. n 

and f' (n-l) = o '»l-p A *[au(n) + bu(n-l)] 
where u(n) is a wh_te noise. 



(A- 22) 



The initial condition is 



. o R - u( 1) 



- 

P x CD 


tt 


1 

t — ■ 

i 


P 2 U) 




0 






— 



4. Dynamic Model of Two-dimensional 2 . . . 

Autocorrelation Function R(m,n) = o R . p R ' in '. p y ' n 



Given the autocorrelation function 



m 



n 



R(m,n) = o R *. p h ,,U| .p v '“' (A- 23) 

the discrete power spectral density is (using eq. A- 5): 



S ( z 1 , z 2 ) = S 1 ^ 2 1 - * S 2^ Z 2^ = 



o R 2 .( 1 -P h 2 K 1 -P v 2 ) 



(1 ‘ p h' I l ) ‘( 1 ' p h‘ 2 l' 1)-(l 'V Z 2)‘U-P V ‘ z 2 ' 1 ) (A " 24) 

1 -1 



Separating S(z^ z 2 ) to H(z^ z 2 )*H(z^ > z 7 ) gives: 

,Tu-P h 2 Hl-P v 2 )' 



H ( z , z 9 ) 



R 



1 i '2 J " TT TT 

(1-Pb-Zj. >(1-P^z 2 X ) 



(A * 2 5 ) 



and the output of such transfer function 



G(z x z 2 ) *■ H ( z ^ ^ z 2 ) *u( z 1> z 2 ) 



y.elds the following recursive equation: 



G ( Z 1 , 2 2^ = -h’ z \ l ' G ( z i > 2 2- ) ' p h* p v‘ z l 1 * z 2 1 * ( ^^ z l, z 2' 



* o v *z: •G(z 1> z 2 ) 



+ a 



R 



i/(l-P h Z Hl-P v 2 )-u(z lf Z 2 ) 

( A- 2 6) 



2 4 8 



The dynamic model in spatial domain will be: 



g(in,n) = P,,*B(m-l,nJ - P,*p v ' g (m- 1 , n- 1) + P v *g(m,n-1) 

+ o R / (1-P h 2 >(1*P V 2 )' ’ u(m,n) C‘ A - 2 7 ) 

for m,n > 1 



Equation (A-27) is a two dimensional first order Markov 
process representation. 

The initial condition will be: 



g(l,l) = o R *u(l^l) 




(A-27) 


g(m,l) = P h * g (it- 1,1) 


+ a R^ 1 ' p h 2 * u ( ra » 1 ) 


m > 1 


g(l.n) = P y • g ( 1 » n- 1) 


+ a R^ 1-P v 2 ‘ud.n) 


n > 1 



where u(m,n) is white noise with variance 1. 

This dynamic model was used by Habibi [Ref. 17] and 
can be characterized in a state space structure [Ref. 18] by: 



r (m+ 1 , n) 


3 


r < i 

_ p hj_ 


• 


* 

r (m,n) 


s (m ,n+l) 

• mt 




0 1 P V 


1 


s (m , n) 

m 



(A- 29) 






0 






T 1 



u (m , n) 



gO,n) =» ( 1 



0 ] • 



r (m,n) 
s (m , n) 



5 . 



Uy nam i c Mud el o 1 Two - d imens i ona 1 
Autocorrelation l : un c t i on R (in , n ) = 

2 | in ! | n | r \ r ' 

O n • D. . D VC O S ( 0.1 • III ) c OS (a) Ml) 

R h v h 7 v v * 



Given the autocorrelation function: 



R(m,n) = a R ~.p h 1 m 1 P v ' n 1 . cos (ui^ *nO • cos (cj v • n) 



(A- 30) 



The discrete power spectral density (using eq. A-10) is: 



S ( Z 1 , Z 2 ^ = S ( z l^ ’ ^ (z 2 ) = 



where 



= H l (z l ) • H^z^ 1 ) * H 2 (z 2 ) • H 2 (z 2 ' 1 ) (A- 31 ) 
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- °R/U-e h 2 Hl-» v 2 )' 
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where 
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The output of the transfer function of eq . (A-32) is: 

G( z. z ? ) = II (.2 2 . ) * U(z, z ? ) 

1,4- 1,4 X,4. 

Calculate G(z. z 7 ) and make the two- diinens ional inverse z 
transform, the dynamic model in spatial domain will be: 

2 

g(m,n) = C v *g(m,n-1) - P v *g(m,n-2) 

+C^g(m-l,n) - C h «C v «g(m-l,n-l) + c h *p y 2 *g(m-l,n-2) 

• D h 2 .g(m-2,n) + P h 2 *C v *g (m-2 ,n-l) - p 2 h • p y ^ g (m- 2 , n- 2) 

+ o r Y (l-p h 2 ).(l-p v 2 ) , .[a h *a v *u(m,n) + a h * a y * u(m, n- 1) 

+ a v *b h «u(m-l,n) + b h -b v *u(m-l,n-l)] (A-34) 

for m,n > 2 

where = 2‘p^‘Cos 

c v = 2 • p v • cos u» v 

and u(m,n) is an uncorrelated random field (white noise). 

The initial condition is: 
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Equation (A-34) is a two dimensional second order Markov 
process representation. 






An example ot an image 100 x 1UU in size, formed by 
equation (A-5a). representing the given autocorrelation func- 
tion eq. (A-30) is shown in Fig. A- 1 In this picture we used 

p. = p = 0.99 and o = l and u> = - 0.143. Comparison of 

h v K n v r 

the correlation 'unction R(m,n) of the simulated image to the 
given image yielded an average error of the correlation factor 
around l ? o . 

The state space structure can be found as follows: 



Let us def ine : 
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By using the following definition: 
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Figure A-2 
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2p v cos • p 4 (z l, z 2 ) * z 2 _1 ' p 3 (z l, z 2 ) 

'll. 

+ z 2 -U (z x z 2 ) 



Converting eq . (A-43) to spatial domain and putting in state 

variables f o nn , we get: 
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where 



and u (m,n) 



g(m,n) = [0 1 0 0] 



is an uncorrelated random 
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APPENDIX B 



BAYES CRITERION 

A. REVIEW OF ONE DIMENSIONAL CRITERION [32] 

In this appendix, a threshold setting algorithm will be 
presented, similar (but not identical) to the Bayes decision 
algorithm. First, the basic decision model in one dimension 
will be reviewed. 

- Given a measurement sample x - 

x = s + n 
s = desired signal 
n = noise 

- There is an unknown parameter 0 related to the measure- 
ment. Example: 

Given a pixel containing noise and "possible" signal, 
as in Fig. B-l(a) - 




Pixel (m) Pixel (m) 

Fig. B- 1(a) 



The "possible signal that has to be detected in pixel (m) is 
shown in Fig. B-l(b) - 



256 



r Si 

Signal = s (m) =< or 
ls 0 

Pixel (m) 

Fig. B-l(b) 

The resulting signal plus noise samples are shown in Fig. 
B-l(c) - 




Measurement samples = x(m) 



i 1 

i i 

I I 

I I 

"T i 



Pixel (m) 



Fig. B-l(c) 



The unknown parameter ("event") 0 - 



0 = < 



S 1 9 1 
or 

s 0 " 9 0 



Given the conditional PDF's 



P[x/0 Q J P[x/0 1 ] 



such as shown in Fig. B-2 - 
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Fig . B- 2 

In order to determine whether the signal in pixel (m) 
is Sq or s^, the following definitions have to be made: 

- Define the following decision errors: 

a 2 Reject when correct 

a - L P[X/0,]dx 
*1 1 

g = Accept when false 

g = L P [X/0] dx 
! 1 

- Define the following costs: 

= Cost of deciding "Sq" when not true 



= Cost of miss 
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Cp^ = Cost of 



deciding 



when not true 



= Cost of false alarm. 

(In the Bayes solution there are two additional cost 
terms which are not required here.) 



Define the risks: 



C 



M 



♦6 



R 



0 



C 



FA 



a 



R 



1 



Given the probability of "event" 0 q =*> ("event" of 
no target) 

Given the probability of "event" 0^ => g^ = l-g Q ("event" 
of target present) 

The total risk can be found by combining all the defini- 
tions above and is given by: 



R = cr • R + cr • R 
°0 0 °1 1 






•0 M 

V. 

°o l m -co 



(B-l ) 



TH 



P[X|0 Q ]dx h- g 



1 FA 



V 



[X | ©J 



d.x 



TH 



Minimizing R with respect to V_,„ ( 
the following threshold equation: 



dR 



W 



TH 



0) yields 



g o C M P[ V TH/0 o ] - gl C fA P[' , TH/0 1 ] (B- 
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lO 



When all the parameters are known, we can solve (B-2) 
for the threshold level. The following two definitions are 
useful : 



. P [X | ©! ] 

P [X | 0 Q j 

A(X) = 



“ A (*) 
likelihood 



ratio 



S o C M 



= k 

threshold ratio 



( B - 3 ) 



( B-4 ) 



Example 

Given a measurement X that can belong either to X Q or X 



where 



X Q = N ( 0, a ) 
§o’ 9 o = 0 



X x = N ( y , a “ ) 



§1 = 9 1 



d 



N is normal PDF. 



Decide whether X belongs to X Q or X-^ 
- Decision method: 



P[X /0] = 



P^/y] = 



/2tto 

1 

/Tr?a 



2 2 
- X / 2a “ 



(X-y) 

2a 2 



k = 



g o C M 

S 1 C FA 



r r a a assumed to be specified. 

u FA’ °o’ 0-1 
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The likelihood ratio is: 



-CX-y) 2 /2o 2 

A(X) = 2 — r~ 

-XV2o Z 



The threshold level is: 

V TH ’ I * T ln K 



= K 

x=v TH 



(B-S) 



The decision algorithm: 

when X > X belong to X^ 

when X < X belong to Xq 

c. Extension to Two Dimensional 
Multivariable Criteria 

Following the one dimensional case, we can 
define the signal to be a vector rather than a scalar. 



The signal is S 
The noise is N 
The measurement is X 



[ Sl s ? •' 


* s I ] 


[n. n ? . , 

± y “ y 


..Hi]. 


[ x 1 , x 2 , ■ ■ 




s T * n t 





For the multivariable normal PDF, the Bayes 
solution requires that the vector X will be filtered first 
by a matched filter and then compared with a given threshold. 

- This conclusion can be demonstrated as follows: 

The multivariable normal PDF is: 
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P (N) 



1 



( B-6) 



(2 it) 1/2 | R 



1 / 2 * Exp [- 2 N n ] 



NN 



P(N) = multivariate normal PDF 



R 



NN 



= covariance matrix of the noise 



= E[N T N] 



(assuming E[N] = 0) 



The noise vector can be presented as follows: 



N T = (X T - S T ) = { 



(X T - S X T ) 



(X T - S n T ) 



or 



A point target can be expressed as : 
S ] _ = [0,0, .. . s (m,n) ,0] = S 
S Q = [0,0, ... 0. .. 0] =0 



The likelihood ratio is given by: 

Exp[4(X - S 1 )R nn ‘ 1 (X - Sl ) T ] 

A (X) — n =j Fn "* 1C 

Hxp[-4(X - - S o rj ( B ‘ 7 ) 

Taking the natural log of ( B- 7 ) will yield the 
required conclusion: 

In k - (X - S)R NN 'hx - S) T ♦ 4x R nn X T 

= X.R _1 .S - C ( B ■ 8 ) 

x r nn b L 1 

where 

1 -IT 

c i * 7 sr nn s 
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Input I 
vector! 



< In k + Ci =o S „ 
x J MatchedI - 






filter 



> lnk + C x ^ S 1 



as described in Fig. B-3. 




1 = target 



0 = no target 



Figure B-3 
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APPENDIX C 



DERIVATION OF ADAPTIVE FILTER ALGORITHM* 



A. ALGORITHM DERIVATION 

At a given pixel, the output of the filter is 



y . = X . T W = W T X. 
J J J 



(C-l) 



Let us select a desirable reference performance, , 
which we seek to achieve at the steady state. For example, 
in the case of point target detection (or estimation) , we 
seek to achieve the best estimate of a point target while 
suppressing the background clutter noise as much as possible 

Let us define the error e. at pixel location j to be 

J y J 



e . = d. -y. = d. - W 1 X. 
J J J J 



(C-2) 



It is the purpose of the adaptive process to adjust the 
weighting coefficients of the linear adaptive filter to 
minimize the mean square error, e . . 

The derivation of the adaptive algorithm is started by 
considering the k-th step of adaptive iteration. The mean 
square error is given by 



£ . 2 = d. 2 - 2d .X . T W, + W, T X. X. T W, 

J J J J k k j j k 



(C- 5) 



Its statistical expected value is obtained by averaging (C-3) 
over an ensemble of identical adaptive filters. 



*This appendix is based on the original derivation of 
B. Widrow [21,22] . 
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E[e. 2 ] . = E[d. 2 ] - 2E[d.X T ]W T + W, T E [X . ♦ X T ] W. 

L 3 J w=wk 3 Jjk k L j jk 

CC-4) 

Let us define a vector P as the crosscorrelation between 
the desired response and the X- vector. 



>T 



= E [ d . X . T ] = E[d.X, ,d.X,,..d.X T ] 

J J J 1 •* J 2-' J I • 

j j J J J J j j 



(C-5) 



the ensemble average 



R = R = E [ X • X T ] = 
xx 1 J 3 



r~ 


XX 
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X x ) (Xp 
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x 1 ) . . . 
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i 



( x i x j ) 

j j 



The mean square error £ ^ can be expressed as 

h - £ [ d / 1 - 2pT 'v \ Tr w k 



C C - 6 ) 



(C-7) 



It is a quadratic function of the coefficients. Adjust- 
ing the coefficients involves descending along a concave 
hyperparaboloidal surface with the objective of reaching its 
minima. Gradient method is commonly used for this purpose. 
The gradient of the mean square error function at W = 
is obtained by differentiating (C-7) 

= - 2 P + 2R ' W, 




w = wk 



(C- 8) 
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The optimal filter W* , generally called the Wiener filter, 



is obtained by setting the gradient to zero and is: 




(C-9 ) 



For some of the subsequent analysis, it is convenient to 
express the mean square error £ k (C-7) and the gradient 
function ( C -8 ) in another form to be derived as follows: 
Substituting ( C-9) into (C-7) will yield: 

'min* E ( d j Z l - W * T P ( c ' 10 ) 

Recombining (C-10) with (C-7) and (C-9 ) will yield: 



^k ^min 



T 

V. R 
k 



V, 



(C-ll) 



where V^- - W* (C-12) 

and the gradient becomes 

V k = 2RV k ( C - 1 3) 

Since R is the correlation matrix, it is positive def- 
ite, and may be expressed in normal form as follows: 

R = Q A Q " 1 ( C - 14) 

where A is a diagonal matrix of The eigenvalues of R. If Q is 

- 1 T 

constructed to be orthonormal, then Q = Q and then ^.^ 5 ^ 
R ' 1 = QA ' 1 Q ' 1 (C-16) 

and the mean square error may be expressed as: 
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(C- 17) 




A new set of coordinates may now be defined as 



v = Q T v = Q" 1 V 



(C-18) 




(C-19) 



and 




(C-20) 



This transformation projects V into prime coordinates. 
Since A is diagonal, the prime coordinate must comprise the 
principal axes of the quadratic mean square error performance 
surface. The gradient expressed in primed coordinates becomes 



B. THE METHOD OF STEEPEST DESCENT 

The objective of the adaptive process is to find a 
solution to (C-9). It should be noted that (C- 9) has exactly 
the same form as the equation in Chapter III for designing 
the spatial MMSE filter. 

In that chapter, nonadaptive filters were developed. 
Solution of (C-9) required a two step process by first 
calculating the statistical properties of the image and of 
the targets, then followed by a solution of a set of linear 
equations . 



V 



k 




(C- 21) 
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In this chapter, adaptive filters are being developed. 

An alternative approach is used, based on the steepest descent 
adaptive method. 

The adaption in this method starts with an initial value 
assigned to the filter -W . The gradient of the mean square 
error is measured and the filter coefficients are modified 
accordingly. This procedure, when repetitively applied, 
causes the error to reduce successively until the filter 
approaches the optimal value. This method can be described 
by the following relation. 

(C" 22) 

where y is a parameter (loop gain) that controls stability 
and convergence rate, is the value of the gradient at a 
point on the error surface corresponding to W = W^. 

An expression for is given by (,015). Substituting 
(013) into (022) will yield 




W 



k+1 



= W, - 2y R V, 



(023) 



Subtracting W* from both sides of (0 23) will yield 



V. 



. + 1 = V k - 2y R V k = (I - 2U R) V } 



( C - 24) 



Equation (C - 2 4) is a linea 
equation whose solution charac 
of the adaptive filter. The s 



r homog 
ter izes 
olution 



enous vector difference 
the dynamic behavior 
to (0 24) i S : 



V k = (I - 2uR) k V Q 



( 0 25 ) 
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It is stable if 





lim (I - 2yR) k = 0 
k-voo 


(C- 26) 


Since 


(.1 - 2y R) = Q Cl - 2yA)Q' 1 


(C-2 7) 


and 


Cl - 2y R) k = Q(I - 2yA) k Q' 1 


(C-28) 


condition 


(C-26) will be satisfied if 






lim (I - 2yA) k = 0 

k->oo 


(.C-29) 


Condition 


(C-29) will be satisfied if 






|1 - 2y X p | < 1 


(C - 5 0 ) 


for p = 


1, 2. . . I 




or - — — > 

A max 


y > 0 


(C-31) 



where ^ is the largest eigenvalue of R 

max ° * xx 



Equation (C-31) gives the stable range for y. 

It can be seen that in primed coordinates, the coefficient 
equation given by 

V' k+1 = (I - 2yA)V’ k (C-52) 

with the solution 

V' k = (I - 2yA) k V o (C- 35) 

For the p-th coordinate we may write 

Vk ‘ C1 ' 2,J> P )k V po (c - 34) 
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Eq. (C-34J represents a geometrical progression 
starting from the initial condition v ' . 

The geometric ratio of the p coordinate is: 



for 




r p f 1 ■ 2v y 



(C- 35) 



Since we prefer to talk usually in terms of 
progression, we can fit an exponential envelope 
with time constant, t . If the unit of time is 
tion cycle, then 



exponential 
to C C - 33) 
one itera- 



r = exp ( 1 /t ) 
P P 




(C- 36) 



In practical applications, y is selected so that t^ >> 1 
and (C-56) can be approximated by the first two terms only. 

Combining (C-36) and (C-35) gives the formula for the p-th 
time constant 



r 



P 



! - 2 y A p 




T P 



1 

2y A p 



( C " 5 7 ) 



Therefore, we can summarize that the transients in the 
coefficients of ’’filter V" are geometric sequences with an 
approximate time constant, r . In the filter W (not in 
prime coordinates) the transients will consist of sum of 
geometric s equenc es , with the total number of time constants 
typically equal to the number of coefficients. 
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Similarly, while the coefficients are varying, the mean 
square error (MSE) varies as well. Substitution of (C-33 ) 
into C C- 20 ) will yield the "learning curve" - the MSE as a 
function of k: 

■ ^jnin * V o' h 1 ' W 2 ** V o IC-38 ) 

As long as condition (C-34 ) and (C-31 ) are met, the adaptive 
process will converge on the minimum point of the MSE surface. 



lim 5 k 
k-^°° 



^min 



( C - 3 9) 



Here again 

sequences 

resulting 



, the relaxation process is a sum of geometric 



where the p-th mode has a geometric ratio of (1 - 2y A ) 

ir 



in a time constant 
_ 1 



fbs< 



4y A . 



= 



(C-42) 



To accomplish the coefficient's adaptation, we have to 
measure or estimate the gradient at each iteration. Reference 
[ 21 ] suggested two methods to obtain the gradient estimate. 
The first method uses differentiation and requires that finite 
perturbations be made on the filter coefficients. The 

second, so-called LMS method, obtains gradient estimates 
directly and without perturbing the nominal coefficients 
adjustment . 

In this work, the second method was selected, based on 
two reasons . 
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1. Practical implementation and utilization of this method 
„ is simpler, compared to the other alternative. 

2. Reference [ 21 ] demonstrated the superiority of the LMS 
algorithm with regard to convergence time and mean square 
error, compared to the other candidate. 



C . THE LMS ADAPTATION ALGORITHM 

The LMS algorithm is an implementation of the steepest 
descent method that employs an efficient gradient estimation 
technique. This method was found applicable to the two di- 
mensional filtering configuration presented earlier in Fig. 



7-3. 



1 . Gradient Estimation, Convergence and Time Constants 
The error, £j , of a filter such as shown in Fig. 7-2 
is given by (C-2). A gradient estimate may be obtained by 
squaring the single value of e- and differentiating it as if 
it were the mean square error: 





r, 2 1 




' " 


/V 


3w^ 




i — 1 

ro 


V. = 
J 


♦ 

* 2 

0 £ j 


= 2tj 


4 




8 w j 




3Wj 




* ^ 




^ m 



npi 3 £ • 

since e- = d. - W A X- and - — 4- = -X.. 

J J j 3w j 



(C-41) 



Substituting (C-41) into (C-42) will yield the LMS algorithm 



W., = W • + 2y£ .X . 

J + 1 J M J J 

where the index k is replaced by j since we 
tation step for each new input vector, i.e. 
data samples. 



(C-42) 

perform an adap- 
a new set of 
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The gradient estimate of (C-41) may be implemented in a 
practical system without further squaring or averaging and 
is simple and elegant. All components of the gradient vector 
are obtained from one input data set X. Since the estimate 
is obtained without averaging, it contains a large component 
of noise. The noise, however, is averaged and attenuated by 
the adaptive process, which acts as a low pass filter in this 
respect . 

It is important to note also that for a fixed value of 
W, the estimate is unbiased, which can be shown as follows: 

E[Vj] = -2E[e.X.] = -2E[d j X j -XjX^W] (C-43) 

From (C-8 ), this expression can be rewritten as 

E[Vj] = -2(P-R4V) = v ( C - 4 4) 

recalling that v is approaching zero at steady state. 

2 . Convergence of the Mean of Filter Coefficients 

Let us now show that the mean value of the filter 
coefficient instead converges to a steady state when the LMS 
algorithm is implemented. Taking the expected value of 
(C-42)j we obtain the following: 

E [W • ] = E[W. ] - 2UE(x.[d -X. T W,]} = 

J +1 J J J J J J 

= [I + 2yR ] .E[W ] - 2y P (C-45) 

In obtaining (C-45), we assumed that W and X are 
uncorrelated. It will be shown later that proper selection 
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of the adaptation loop gain y will cause the filter coeffi- 
cients to vary very slowly compared with the value of the 
noise (X) correlation coefficients. 

E [x T w] - E [x T ] -W 

where W is considered almost as a constant relative to 

'j' T T 

E[x ]. Under these assumptions, E[x *w] - E[x ]*E[w], i.e., 

X and W can be considered uncorrelated. 

Referring back to equation (C-45), with an initial 

filter coefficient, W , j + 1 iterations of equation (C-45) 

yield 

E[w j + 1 ] = [I + 2yR xx ] j + 1 W o - 2y J [I + (C-46) 

Eq. (C-46) may be put in diagonal form by using the normal 
expansion of R , that is, 

AA 

R xx = Q’ 1 

The eigenvalues are all positive, since R is positive 
definite. Equation (C-46) may now be expressed as 

E[w. + 1 ] = [I + 2yQ _1 AQ] ] + 1 W - 2u.Z [1 + 2 Q - 1 A Q] i P = 

J ± u 1=0 

= Q" 1 [I + 2yA] j + 1 Q-W o - 2yQ" 1 J o [I + 2y A] 1 Q-P (C . 47o 

Consider the diagonal matrix [I + 2y A ] . As long as its diagonal 
terms are all of magnitude less than unity 

lim [ I + 2y A] +1 -> 0 

i -boo 
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and the first term of (C-47) vanishes as j-*» . 

The second term in (C-47) converges to a non-zero 

limit. The summation factor - 

Z [I + 2yA] i 

i=o 

j i 1-1 

becomes lim z [I + 2yA] = - y- A 

j-»-oo i=o ^ 

where the formula for the sum of the geometric series has 
been used, that is: 



ii 0 (1 + 2|jA p )1 ■ 1-cWp) = 2 U A p 



Thus, on the limit, equation (7-49) becomes 



lim E [w . , ] = Q _ 1 AQ-P = R -P 

■ -r-oo J ^ 1 AA 



1 

Comparison of this result with (C-9.) shows that as the 
number of iterations increases without limit, the expected 
value of the filter coefficient converges to the coefficients 
of the Wiener filter. 

The convergence is assured if, and only if, y is set 
within certain bounds. Since the diagonal terms of [I + 2yA] 
must all have magnitude less than unity, and since all eigen- 
values in A are positive, the bounds, on y, are given by 

I P A max i 



'•max 



< y < 0 
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where Amax is the maximum eigenvalue of R . This convergence 

XX 

condition can be related to the total input signal power as 
follows : 



Amax < trace R 



xx 



I 



(C-48) 



Therefore, the bounds 
-1 



trace [R ] - E[xJx.] = £ E[x.^] = total input 

XX J 3 i=l 1 power. 

for satisfactory convergence are 



< y < 0 . 



2 E[x. ] 
i= 1 1 



In practical applications, when precise adaptation is required, 
y is usually selected such that 

(C -49) 




3 . The Noise of Gradient Estimation 

Let it be assumed that the adaptive process, using a 
small loop gain y, has converged to a steady state near the 
minimum of the MSE (5). The gradient estimation noise of the 
LMS algorithm at the minimum, where the true gradient is zero, 



is the 


gradient 


es t imate 




N. = 


/>. 

V- = -2c 




3 


3 


since 


II 

< o 


<1 

+ 




3 


3 3 


where 


V. = 


the true 




3 






II 

' t > 


gradient 



J 3 



(C - 5 0 ) 
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The covariance of this noise is given by 



cov[N-j = E[N.N.^] = 4E [ £ - 2 x • x .^ ] 
L J L 3 J 1 L 3 3 3 



(C-51) 



It is known from Wiener filter theory that when Wj = W* [the 
optimal filter] the error £. is uncorrelated with the input 
vector x.. If one assumes that £. and x. are Gaussian, they 
will be independent, as well, at the minimum. Under these 
conditions (C-51) becomes 



cov [N . ] = 4E [ e • 2 ] * E [ X . X - T ] = 4C . R 
L j J 1 j J L j j J mm xx 



(C - 52) 



In primed coordinates, the covariance is 

cov[N'.] = Q ^cov[N.]Q = 4£ . A 
L j x J mm 

4. Noise in Filter Coefficients 



(C- 53) 



Adaptation based on noisy gradient estimates results 
in noise in the filter coefficients. The method of steepest 
descent with ideal gradients is represented by (C-29). With 
estimated gradients, this equation may be rewritten as: 

V • , -j = V. + U (-V.) = V. + (-V. - N.) (C- 54) 

J + 1 J ‘ ^ J 3 ^ 3 J K J 

~ A 

where the estimated gradient is V. = V. + N-. 

J J J 

Substituting (C-13) and combining terms yields 



V 



j + i - Ci • 2uR xx )v j 



(C - 5 5 ) 



This is a first order difference equation with a stochastic 
driving function of -yN. . Projection into primed coordinates 
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may be accomplished by premultiplying both sides of (C-S5) 

by Q" 1 . 

V' j + 1 = (1 - 2yA)V' ;j - y N j (C- 56 ) 



In steady state, after initial adaptive transients have died 
out, Vj' undergoes a stationary random process in response to 
the stationary driving function, -yNA . Since there is no 
cross coupling between terms and the components of N4 are 
mutually uncorrelated, the components of V4 will also be 
mutually uncorrelated, and the covariance matrix of N'. will 
be diagonal. To find the matrix, one first multiplies both 
sides of (c-56 ) by their transposes. 

V' j + 1 V. + i T = (I - 2yA)VyV.' T (I - 2yA) +p 2 N' j N j ' T 

- y Cl - 2yA)V j 'N'. T - y N j , V j ,T (I - 2yA) 

(C-S'7) 



Taking expected values of both sides yields 



cov[Vb] = (I - 2yA)cov[VV] (I - 2yA) + y 2 cov[N^.'] 

( C-58 ) 

(Note that VA ' Nj 1 are uncorrelated since VA is affected only 
by gradient noise from previous adaptive cycles.) 

Combining terms further yields: 

cov [V^ ' ] = y 2 (4yA - 4y 2 A 2 ) _1 cov [N ' j ] ( C-59) 

In practical applications, when y is small, 
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y A << I 



Neglecting the squared terms in (C-S9) 



cov [Vj '] = ^ A" x cov[N9 ] 



Combining (C-59 ) with (C-52), one may write 



cov [V h ] = 4 T A A C4 ?min A) = y? min I 



-1 



(C-60 ) 



The covariance of Vj can be obtained from (C-60) by using 

(C-16 ) and (c-18 ) • 

cov [V . ] = E[QV.' v.' T Q" 1 ] = Q-y£ . IQ" 1 = y£ • ‘I 

l jj j j x j x x H ^min 

C C * 61 ) 

5 . Misad j us tment 

Without noise in the weighting vector, adaptation by 

the method of steepest descent would converge to a steady- 

state solution at the minimum point of the MSE surface. The 

MSE and the minimum will be £ min - Noise on the weighting 

vector, however, tends to cause the steady-state solution to 

vary randomly about the minimum. The result is an excess mean 

square error, a mean square that is greater than 

An expression for mean square error in terms of V' is 

given by (C-20), where the excess mean square error is 
T 

V. ’ AV.' , The average excess mean square error is 
JJ 

E[V ' T AV '] = E X p E[v - rj - P5 min E Xp 

J J P=1 P=1 



where, according to (C-61), 
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E[(v ' ) 2 ] = U5 min for all P 

A useful parameter in the design of adaptive processes 
is the misad j us tment M, which is defined as the average excess 



mean square error divided by the minimum mean square error. 

, T 

“ v J \ 

l C - 6 2) 



M 



A E[V'. T AV.]] 



■mm 



The misadjustment is a dimensionless measure of the difference 
between adaptive performance and optimal Wiener performance 
as a result of gradient noise. The misadjustment is thus 
given by 




(C - 65) 

This useful formula may be re-expressed in a manner that 
allows one to perceive the relationship between misadjustment 



and rate of adaptation. According to (7-42) one may write 

1 ' ' ' (C-64) 



1 



» tr R xx = T Ji, — 

P = 1 D 



= zCc4-) 



ms e 



p av 
^mse 



The misadjustment may thus be written 

I 
4 



M - t ( ) 



( C - 6 5) 



' p av 
^mse 



where I is the number of filter coefficients and x is 

Pms e 

the MSE time constant, related to the learning curve. When the 
eigenvalues are equal, 

( C - 6 6) 



M = ^(^-) 



mse 
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Since trace is the total power of the input vector X, 

which is generally known, one can use ( C - 6 3 ) in choosing a 
value of that will produce a desired value of M. Combin- 
ing C C -66) with CC -63) will yield a general formula for the 
time constant of the learning curve with equal eigenvalues. 

x = I (C- 67) 

mse 4u tr R YY 

AA 

This formula is also a good approximation in many cases when 
the eigenvalues of R YY are unequal. 
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